Matrix Computations Second Edition Gene H. Golub Charles F. Van Loan 1989 Голуб Дж., Ван Лоун Ч. Матричные вычисления: Пер. с англ. - Мир, 1999. - 548 с., ил. Книга известных американских математиков-вычислителей представляет собой удачное сочетание учебного пособия и справочника по методам численной алгебры. Изложение сжатое, в рецептурной форме, без доказательств. Книгу отличают методические достоинства: каждый раздел содержит задачи для читателей-студентов и обзор научной литературы - для специалистов. Оглавление Предисловие редактора перевода 5 Предисловие к первому изданию 7 Предисловие ко второму изданию 11 Глава 1. Умножение матриц 16 Глава 2. Матричный анализ 56 Глава 3. Линейные системы общего вида 87 Глава 4. Линейные системы специального вида 127 Глава 5. Ортогонализация и метод наименьших квадратов 181 Глава 6. Параллельные матричные вычисления Глава 7. Несимметричная проблема собственных значений 299 Глава 8. Сииметрична проблема проблема собственных значений 368 Глава 9. Методы Ланцоша 426 Глава 10. Итерационные методы для линейных систем 452 Глава 11. Функции от матриц 482 Глава 12. Специальные разделы 501 Предметный указатель 536 Глава 1. Умножение матриц 16 1.1 Основные алгоритмы и обозначения 16 1.2 Учет структуры матрицы 29 1.3 Блочные матрицы и алгоритмы 36 1.4 Некоторые аспекты векторно-ковейерных вычислений 45 Глава 2. Матричный анализ 56 2.1 Основные сведения из линейной алгебры 56 2.2 Векторыне нормы 59 2.3 Матричные нормы 61 2.4 Матричные вычилсления с конечной точностью 65 2.5 Ортогональность и сингулярное разложение 73 2.6 Проекции и CS-разложение 77 2.7 Чувствительность квадратных систем к возмущениям 80 Глава 3. Линейные системы общего вида 87 3.1 Треугольные системы 87 3.2 LU-разложение 92 3.3 Анализ ошибок округления в методе исключения Гаусса 102 3.4 Выбор ведущего элемента 106 3.5 Уточнение и оценивание точности 119 Глава 4. Линейные системы специального вида 127 4.1 Разложение вида LDM LDL 127 4.2 Положительно определенные системы 132 4.3 Ленточные системы 141 4.4 Симметричные неопределенные системы 150 4.5 Блочные трехдиагональные системы 159 4.6 Системы Вандермонда 166 4.7 Теплицевы системы 171 Глава 5. Ортогонализация и метод наименьших квадратов 181 5.1 Матрицы Хаусхолдера и Гивенса 181 5.2 QR-разложение 195 5.3 Задача наименьших квадратов: случай полного ранга 205 5.4 Другие ортогональные разложения 215 5.5 Задача LS неполного ранга 221 5.6 Взвешивание и итерационное уточнение 230 5.7 Квадратные и недоопределенные системы 234 Глава 6. Параллельные матричные вычисления 6.1 Операции на распределенной памяти 238 6.2 Операции на общей памяти 252 6.3 Параллельное умножение матриц 262 6.4 Кольцевые процедуры разложения 273 6.5 Сеточные процедуры разложения 281 6.6 Методы разложения на общей памяти 290 Глава 7. Несимметричная проблема собственных значений 299 7.1 Свойства и разложения 300 7.2 Теория возмущения 307 7.3 Степенные итерации 316 7.4 Хессенбергова форма и вещественная форма Шура 324 7.5 Практический QR-алгоритм 334 7.6 Методы вычисления инвариантных подпространств 343 7.7 QZ-метод для Ax = lyambdaBx 353 Глава 8. Сииметрична проблема проблема собственных значений 368 8.1 Математические основы 368 8.2 Симметричный QR-алгоритм 376 8.3 Вычисление SVD 383 8.4 Некоторые специальные методы 392 8.5 Методы Якоби 399 8.6 Метод разделяй и властвуй 412 8.7 Более общие проблемы собственных значений 418 Глава 9. Методы Ланцоша 426 9.1 Выводы свойства сходимости 426 9.2 Практические процедуры Ланцоша 433 9.3 Приложения и обобщения 442 Глава 10. Итерационные методы для линейных систем 452 10.1 Стандартные итерации 452 10.2 Методы сопряженных градиентов 461 10.3 Сопряженные градиенты с предобусловливанием 471 Глава 11. Функции от матриц 482 11.1 Спектральные методы 482 11.2 Аппроксимационные методы 488 11.3 Матричная экспонента 495 Глава 12. Специальные разделы 501 12.1 Задача наименьших квадратов с ограничениями 501 12.2 Выбор подмножеств при помощи SVD 509 12.3 Общая задача наименьших квадратов 514 12.4 Сравнение подпространств при помощи SVD 518 12.5 Модифицированные задачи на собственный значения 523 12.6 Модификация QR-разложения 528 Предметный указатель 536 Предисловие редактора перевода Вниманию читателя предлагается новая книга по линейной алгебре. И, естественно, возникают вопросы: "Что же действиетльно новое можно найти в ней и чес конеретно она отличается от других книг по линейной алгебре?" Как никакая другая ветвь математики, ленийная алегебра самым тесным образом переплелась с многочисленными приложенеиями, являясь либо предметом, либо инструментом исследований. С основами линейной алгебры с той или иной мере знаком каждый человек, соприкасающийся с вычислениями и тем более имеющий дело с решением больших задач на современной суперЭВМ. Круг лиц, интересующихся линейной алегброй необычайно широк. Он весьма неоднороден по уровню квалификации и включает в себя студентов, аспирантов, молодых специалистов, а также специалистов высшей квалификации. Все эти обстоятельства делают очень трудным написание книг по линейной алгебре. Желание отразить последние достижения и сделать при этот изложение достаточно простым - всегда противоречивы. Член Национальной Академии наук США, профессор Стэнфордского университета Дж. Голуб и профессор Корнельского университета Ч.Ван Лоан являются специалистами в области линейной алгебры, известными не только с США, но и во всем мире. Их имена хорошо знакомы и специалистам в нашей стране. Несколько лет назад они предприняли попытку описать в доступной форму последние достижения в области численных методов линейной алгебры. Эта попытка воплотилась в книгу "Матричные вычисления", первое издание которой было встречено с интересом. За последние пять лет авторы провели значительную работу по улучшению содержания книги: в частности, были улучшены многие доказательства, включен материал, касающийся параллельных вычислений. Результатом этой работы было второе издание книги "Матричные вычисления". Именно оно и предлагается вниманию читателя. Авторам удалось успешно сочетать строгость и простоту изложения. Книга читается довольно легко. в ней имеется много нового материала, который еще не публиковался в систематизированном виде на русском языке. Поэтому она будет полезна практически всем специалистам в области вычислительной математики: студентам и аспирантам - как пособие для освоения основ линейной алгебры, прикладникам - как источник рецептов решения задач, опытным специалистам - - как пища для размышлений. Последнее обстоятельство особенно важно для нашего читателя, так как представленный в книге материал хорошо отражает уровень исследований за рубежом в области численных методов линейной алгебры. Объем второго издания книги "Матричные вычисления" увеличен более чем на треть по сравнению с первым. Введение нового материала во многом объясняется быстрым ростом исследований в области параллельных вычислений. Литература по этому направлению плохо систематизирована. Нет еще установившихся методов исследования. Не всегда даже ясно, что в параллельных вычислениях является главным, а что - второстпенным. Поэтому попытка авторов навести определенный порядок очень важна. Наличие систематизированных обозначений, методов исследования и источников информации позволяет лучше понять сильные и слабые стороны предметной области. Для всех, кто занимается параллельными вычислениями, сейчас это очень нужно. Для удобства читателя списки литературы помещены по главам. Добавим к сказанному, что в книге имеется много интересных решений и находок в отдельных доказательствах. Все это можно увидеть, познакомившись с текстом поближе. Хочется надеятся, что в ряду многих книг в области численных методов линейной алгебры предлагаемая книга Дж. Голуба и Ч. Ван Лоана займет достойное место. В заключение отметим, что перевод книги выполнен группой переводчиков: гл. 7, 8, 11 перевел Ю.М.Нечепуренко, гл. 1, 2 - А.Ю.Романов, гл. 3,4,12 - А.В.Собянин и, наконец, гл. 6, 9 и 10 - Е.Е. Тыртышников. В.В.Воеводин Предисловие к первому изданию Мы придерживаемся той точки зрения, что предназначение численного анализа - обеспечить научную общественность эффективными библиотеками программ. Эта деятельеность особенно интересна тем, что ее участникам необходимо знание знание как математики, так и вычислительной техники. В самом деле, разработка программ высокого качества требует и понимания математического существа решаемой задачи, и чутья к конструированию алгоритмов, и осознания специфики вычислений с конечной точностью. Цель этой книги - снабдить читателя нужными навыками в той мере, в которой они относятся к матричный вычислениям. Начиная с середины 50-х годов в этой области был достигнут огромный прогресс. Это подтверждается существованием высококачественныъх программ для решения систем линейных уравнений, многих задач, решаемых методом наименьших квадратов, и задачи на собственнце значения. Типичный пример - подпрограммы пакетов Eispack и Linpack, широкое использование которых привело к повышению уровня разработки алгоритмов во многих прикладных областях. Используя модули из этих пакетов в качетве строительного материала, ученые и инженеры могут сооружать более сложные, специально настроенные на их нужды, программыне средства. Такое положение вещей стимулирует написание хорошо стркуктурированных программ - благоприятная тенденция в условиях, когда все больше научных результатов формируется в виде программ. Влияние численной линейной алгебры ощущается и в другом. Наши добрые привычки - склонность полагаться на ортогональные матрицы, осознание роли чувствиетльности задачи к возмущениям, тщательное рассмотрение ошибок округления - были восприняты во многих других областях исследований. Хороший пример - рост использования сингулярного разложения (SVD) как аналитического инструмента многими статистиками и инженерами в области управления. Работающие в этих областях специалисты переформулируют многие теоретические концепции на языке SVD, и в результате им становится гораздо легче реализовывать свои идеи в условиях влияния ошибок округления и неточных данных. Дальнейшие свидетельства растущего влияния численной линейной алгебры можно обнаружить в области разработки аппаратных средств ЭВМ. Последние продвижения в области арифметики с плавающей точкой и параллельных вычислений. Мы написали эту книгу, чтобы охватить эту очень интересную и быстро расширяющуюся область с единых позиций. Многое было сделано после публикации в 1965 г. монуметального трактата Уилкинсона "Алгебраическая проблема собственных значений". Многие из этих современных достижений уже отражены в обзорных статьях и в известных специальных монографиях, таких как "Численное решение задачи наименьших квадратов" Лоусона и Хенсона, "Симметричная проблема собственных значений" Парлетта. Мы считаем, что пришло время для синтеза всего этого материала. В этом отношении мы рассматриваем "Матричные вычисления" как всеобъемлющую, несколько более продвинутую версию книги Стьюарта "Введение в матричные вычисления". Мы расчитываем на три категории читателей: студенты-выпускники технических специальностей, ученые и инженеры-вычислители, наконец, наши коллеги по численному анализу. Для каждой из этих групп у нас в книге кое-что припасено. Для студентов (и их наставников) мы включили занчительное количество упражнений. Многие из них носят вычислительный характер и могут стать ядром программного проекта. Наш собственый опыт преподования на основе этой книги говорить о том, что задания, использующие Eispack и Linpack, существенно оживляют ее содержание. Также успешно использовалось совместно с этим текстом система Matlab - простая в обращении и выполняющая вычисления с матрицами. Для инженеров и ученых, желающих в процессе своей работы пользоваться этой книгой как справвочником, мы постарались свести к минимуму взаимозависимочть отдельных глав. Кроме того, в конце практический каждого раздела мы приводим аннотированную библиографию, чтобы ускорить поиск дополнительного материала по любому заданному вопросу. Для наших коллег по численному анализу мы щедро сопровождаем наши описания алгоритмов элементами теории возмущений и анализом ошибок. Мы хотим снабдить эту категорию читателей достаточным количеством подробностей, чтобы они могли ставить и решать матиричные задачи, возникающие в их собственной деятельности. Часто толчком к исследованиям в области численной линейной алгебры служат работы, проводимые в других областях численного аналиоза. Так, например, некоторые из лучших методов решения разреженных линейных систем были разработаны исследователями, занимавшимися численным решением уравнений в частных производных. Подобным же образом развитие приемов модификации различных матричных разложений было стимулировано работами по квазиньютоновским методам. Эта книга писалась шесть лет. Название менялось несколько раз: (1) "Завершающий курс матричных вычислений", (2) "Прикладные матричные вычисления", (3) "Матричные вычисления: специальный курс". Первое название, помимо некоторой претенциозности, вводит в заблуждение. Мы не претендуем на исчерпывающее изложение матричных вычислений. В частности, мы искючили многие темы из животрепещущей области вычислений с разреженными матрицами просто потому, что не хотели погружаться в теорию графов и структур данных. Второе название не подходит, потому что мы не останавливаемя в сколько-нибудь значительной степени на приложениях. По большей части рассматирваемые матричные задачи мы принимаем как данность - их происхождение не прослеживается. Мы согласны с тем, что в педагогической точки зрения это является недостатком, однако опытный преподователь сумеет его компенсировать. А потом, мы подозреваем, что многие из наших читателей и сами будут обладать достаточным опытом, так что чрезмерная мотивировка им и не нужна. Наконец, с послденим названием мы расстались ввиду того, что книга ве же содержит вводный материал. Мы решили включить элементарные разделы и ради полноты, и потому, что мы думаем, что наш подход к основам предмета заинтересует преподавателей начальных курсов численного анализа. О чем же тогда эта книга, если она нполна, не очень-то прикладная и не слишком узко специальная? Ответом на этот вопрос должен послужить краткий обзор содержания. В первых трех главах закладываеются необходимые основы. Приводится обзор матричной алгеборы, фиксируются некоторые ключевые алгоритмы. Темп изложения достаточно высок. Читатели, у которых возникают трудности при решении задач из этих первых глав, безусловно, столкнутся с ними и в остальной части книги. Значительная часть текущих исследований в численной линенйной алегебре сосредоточена на задачах, в которых матрицы имеют специальную структуру, скажем являются большими и разреженными. Искусство извлечения выгоды из наличия структуры - центральная тема гл. гл. 5, где описываются многие методы рашения систем линейных уравнений специального вида. В гл. 6 подхвачено другое течение - возросшее стремление полагаться на ортогональные матрицы. Мы обсуждаем несколько методов ортогонализации и показываем, как применить их в методае наименьших квадратов. Особое внимание уделяется случаю неполного ранга. Центральное место гл. 7 занимает всегмогущий QR-алгоритм для решения несимметричной проблемы собственных занчений. Наш методичный вывод должен помочь снять налет таинственности с этого важного метода. Мы останавливаемся также на вычислении инвариантных подпространств и обобщенной проблеме собственных значений. В гл. 8 мы продолжаем обсуждать проблему собственных значений, сосредоточиваясь на важном симметричном случае. Сначала мы описываем симметричный QR-алгоритм, а затем показываем, как симметрия позваоляет построить несколько альтернативных вычислительных процедур. До самого этого места в книге гаше рассмотрение разреженности носит отрывочный характер. В гл. 5 обсуждается решение ленточных линейных систем, в гл. 7 описан метод одновременных итераций, в гл. 8 - итераций Рэлея и т.д. Главы 9 и 10, однако, польностью посвящены решению задач с разреженными матрицами. Темами обсуждения являются метод Ланцоша и родственный ему метод сопряженных градиентов. Мы показываем, как эти важные алгоритмы могут быть использованы для решения многих разреженных задач на собственные значения, по методу наименьших квадратов, систем линеных уравнений. Цель двух послдених глав - проиллюстрировать широкую применимость представленных в книге алгоритмов. В гл. 11 речь идет о задаче вычисления функции от матрицы, что часто требуется делать в приложениях теории управления. Глава 12 содержить подборку матричных задач, из которых несколько демонстрируют мощь сингулярного разложения. Практическая и теоретическая ценность этого разложения является, пожалуй, сквозной темой в этой книге. Действительно, его алгомитрические и математические свойства играют ключевую роль почти в каждой главе. Во многих отношениях данную книгу можно рассматривать как красочное развитие манускрипта нашего коллеги Алена Клайна: "Все, Что Вы Хотели Знать о Сингулярном Разложении (но боялись спросить)". Пора сказать несколько слов о ссылках на имющиеся программы. Мы ориентируемся на Eispack и Linpack, и почти каждая подпрограмма из этих пакетов так или иначе упоминается в тексте. Кроме того, многие технические доклады, ссылки на которые мы даем в аннотированных библиографиях, являются на самом деле описаниями программ. Следует подчеркнуть, однако, что с большей частью этих программ мы напрямую не работали (исключая Eispack и Linpack), а посему - caveat emptor (лат.)). Предисловие ко второму изданию Две причины побудили нас написать переработанную и дополненную версию "Матричных вычислений". Во-первых, после пяти лет преподавания по прежнему изданию и многочисленных замечаний наших коллег стало ясно, что многое из написанного можно было написать лучше. Во иногих доказательствах и выводах была достигнута большая ясность. Кроме того, мы приняли стилизованный вариант системы Matlab, облегчающий "матрично-векторное мышление". Другая новая черта второго издания - доведение рубрикации материала до уровня пдразделов. Это должно облегичить пользование книгой как учебником и как справочником. Вторая причина появления нового издания связана с расцветом параллельных матричных вычислений. Многопроцессорные ЭВМ революционизируют роль вычислений в науке и технике, и важно, чтобы вклад численной линейной алгебры в этот процесс было документально оформлен. Мы это делаем на машинно-назависимом уровне, подчеркивая общие адгоритмические идеи в противовес конкретным реализациям. Как мы говорим в новой главе о параллелльных матричных вычислениях, эта область очень подвижна, и литература наводнена исследованиями отдельных случаев. Несмотря на это, в области разработки алгоритмов можно набрать дюжинуновых идей, играющих ключевую роль, которые, по всей вероятности, пригодятся нам еще в течение длительного периода и созрели для обсуждения на уровне учебника. Вот обзор новшеств второго издания: В гл. 1 (Умножение матриц) мы вводим обозначения и фундаментальные понятия на примере умножения матриц. Большое внимание уделяется приемам работы с блочными матрицами; добавлен также новый раздел о векторно-конвейрных вычислениях. В гл. 2 (Матричный анализ) мы собрали математические основы, необходимые для вывода и анализа алгоритмов рашения линейных систем и задачи наименьших квадратов. В гл. 3 (Линейные системы общего вида) и 4 (Линейные системы специального вида) были дополнены новым материалом о высокопроизводительной организации гауссова исключения и разложения Холецкого. Акцент делается на реализации, богатые матрично-векторными и матрично-матричными умножениями. Был также добавлен подраздел о положительно полуопределенных матирцах. В гл. 5 (Ортогонализация и наименьшие квадраты) появлися новый материал о блочных преобразованиях Хаусхолдера. Мы также изменили порядок изложения с целью отделить обсуждение ортогональных разложений от решения задачи наименьших квадаратов. В гл. 6 (Параллельные матричные вычисления) мы, пспользуя операцию gaxpy, вводим понятия вычислений над распределенной и общей памятью. Для этих двух моделей параллельных вычислений мы обсудаем организацию перемножения матриц и различных матричных разложений. В гл. 7 (Несимметричная проблема собственных значений) мы добавили новый подраздел о приведении к блочно-хессенбергову виду. В гл. 8 (Симметричная проблема собственных значений) внесены два более фундаментальных изменения. Раздел о методе Якоби был полностью переработан в свете последних достижений в параллельных вычилениях. Был также добавлен раздел о новом высокопараллельном алгоритме типа "разделяй и властвуй" для трехдиагональных матриц. В гл. 9 (Метод Ланцоша) мы добавили подраздел о методе Арнольди. В гл. 10 (Итерационные методы решения линейных систем) расширено обсуждение метода сопряженных градиентов с переобуславливанием. Единственное существенное изменение в гл. 11 (Функции от матриц) и гл. 12 (Специальные вопросы) - это новый подраздел в $ 12.6 о гиперболических модификациях. О сокращениях Следующие книги часто цитируются в нашем тексте: SLE Forsythe G.E. and Moler C. (1967). Computer Solution of Linear Algebraic Systems. [Имеется русский перевод: Форсайт Дж., Молер К. Численное решение систем линейных алгебраических уравнений. 1696 ] SLS Lawson C.L. and Hanson R.J. (1974). Solving Least Squares Problems. [Имеется русский перевод: Лоусон Ч., Хенсон Р. Численное решение задач метода наименьших квадратов. 1986.] SEP Parlett B.N. (1980). The Symmetric Eigenvalue Problem. [Имеется русский перевод: Парлетт Б. Симметричная проблема собственных значений. Численные методы. 1983.] IMC Stewart G.W.(1973). Introduction to Matrix Computation. AEP Wilkinson J.H.(1965). The Algebraic Eigenvalue Problem. [Имеется русский перевод: Уилкинсон Дж. Алгебраическая проблема собственных значений. 1970.] При ссылках на пакеты программ Linpack и Eispack имеются в виду соотвествующие описания: Smith B.T., Boyle J.M., Ikebe Y., Klema V.C., and Moler C.B.(1970). Matrix Eigensystem Routine: EISPACK Guide, 2nd ed, Smith B.T., Boyle J.M., Dongarra J.J., and Moler C.B.(1972). Matrix Eigensystem Routine: EISPACK Guide Extension. Dongarra J., Bunch J.R., Moler C.B., and Stewart G.W.(1978). LINPACK User Guide. Алгольные версии многих подпрограмм из пакетов Linpack и Eispack собраны в книге HACLA Wilkinsom J.H. and Reinsch C., eds.(1971). Handbook for Automatic Computation. Vol. 2, Linear Algebra. [Имеется русский перевод: Уилкинсон Дж., Райнш К. Справочник алгоритмов на языке АЛГОЛ. Линейная алебра. 1976.] Во время написания данной книги (1989), программы, реализующие многие из представленных в ней алгоритмов, можно получить, отправив электронную почту по любому из адресов: netlib@anl-mcs.arpa netlib@research.att.com research!netlib Типичные запросы: send index send index for linpack send svd from eispack Этот сервис электронной рассылки описан в Dongarra J.J. and Grosse E.(1978). "Distribution of Mathematical Software Via Electronic Mail", Comm. ACM 30, 403-407 Сообщим также, что основная библиография этой книги (в формате LATEX) может быть получена по netlib. В каждом случае мы настоятельно рекомендуем включение вычислительных заданий, предполагающих использование высококласных программных пакетов, таких как Linpack и Eispack. Материал может быть еще больше оживлен использованием Matlab - простой в обращении системы, выполняющей вычисления с матрицами. В этой связи мы рекомендуем сопровождающее ее руководство Coleman T.F. and Van Loan C.F.(1988). Handbook for Matrix Computations. Оно содержит главы по Фортрану-77, по элементарным подпрограммам линейной алгебры (BLAS), по Linpack и Matlab. Наконец упомянем новый учебник Hager W.W.(1988). Applied Numerical Linear Algebra. который, возможно, больше, чем "Матричные вычисления", подойдет для начинающих студентов. Глава 1 Умножение матриц Систематическое изучение матричных вычислений начинается с задачи умножения матриц. Хотя математически эта задача очень проста, она таит немало богатств с вычислительной точки зрения. Мы начинаем с рассмотрения в $ 1.1 нескольких возможных способов организации матричного умножения. Заблаговременно вводится язык блочных разбиений матрицы, который мы используем, чтобы охарактеризовать несколько линейно-адгебраических уровней вычислений. Если матрицы обладает какой-то структурой, то обычно удается использовать это обстоятельство. Наприер, для хранения симметричной матрицы общего вида, умножение на вектор матрицы, в которой много нулевых элементов, может потребовать значительно меньше времени, чем вычисление произведения матрицы общего вида и вектора. Эти вопросы обсуждаются в 1.2. В 1.3. вводятся обозначинея для блочных матриц. Блочная матрица - это такая матрица, элементами которой также являются матрицы. Это очень важное понятие как с теоретической, так и с практической точки зрения. С теоретической стороны блочные обозначения позволяют получить очень сжатые доказательства для вожнейших матричных разложений, являющихся краеугольным камнем численной линейной алебры. С вычислительной точки зрения блочные алгоритмы важны, поскольку в них много матричных умножений, а эта операция особенно хорошо реализуется на многих высокопроизводительных ЭВМ новых архитекрур. Эти новые архитектуры требуют, чтобы разработчик алгоритмов уделял обменам с памятью не меньше внимания, чем количеству арифметических действий, что вносит в научные вычисления целое новое измерение. Мы иллюстрируем это в 1.4., где мы рассматриваем такие важнейшие аспекты векторно-конвейерных вычислений, как длина вектора, шаг выборки, количество векторных загрузок/записей, а также степень повторного использования вектора. 1.1. Основные алгоритмы и обозначения Матричные вычисления строятся на иерархии операций линейной алгебры. Скалярные произведения состоят из скалярных операций сложения и умножения, умножения матрицы на вектор составлено из скалярных произведений, перемножение матриц сводится к набору умножений матрицы на вектор. Все эти операции могут быть описаны в алгоритмическом виде или на языке линейной алгебры. Наша главная цель в этом разделе - продемонстрировать, как эти два способа дополняют друг друга. Попутно мы вводим обозначения и знакомим читателя с тем стилем мышления, на котором основываются методы матричных вычислений. Дискуссия вращается вокруг премножения матриц, которое, как мы показываем, может быть организовано несеолькими способами. 1.1.1. Обозначения для матриц Пусть R - множестсо вещественных чисел. Мы обозначаем через R(m x n) векторное просранство всех вещественных (m x n)-матриц: a(11)......a(1n) . . А = a(ij) = . . . . а(m1).......a(mn) 1.1.2. Операции над матрицами Основные манипуляции с матрицами включают: транспонирование С = A(T) => c(ij) = a(ji) сложение C = A + B => c(ij) = a(ij) + b(ij) умножение матрицы на число C = aA => c(ij) = a x a(ij) умножение матрицы на матрицу R(m x r) x R(r x n) -> R(m x n) C = AB => c(ij) = sum(k=1,r)( a(i,k) b(k,j) ) Это те кирпичики, из которых строятся матричные вычисления. 1.1.3. Обозначения для векторов Заметим, что мы отождествляем R(n) с R(n x 1), так что элементы R(n) - это векторы столбцы. С другой стороны R(1 x n) - состоит из вектор-строк. Если x - вектор-столбец, то y = x(T) - вектор-строка. Несколько "подозрительно" выглядит операция над векторами, называемая внешним произведением: C = xy(T), x -> R(m), y -> R(n) Но это совершенно законное матричное умножение, поскольку количество столбцов в x совпадает с количеством строк в y(T). То есть число столбцов вектор-столбца равно 1, как и количество строк в вектор-строке равно 1. Если C = xy(T), то c(ij) = x(i)y(j), так что, например, [1] 4 5 [2] [4 5] = 8 10 [3] 12 15 1.1.4. Операции над векторами Существуют четыре основные операции над векторами. Если a -> R, x -> R(n) y -> Z(n), то мы имеем умножение вектора на число z = ax ( z(i) = ax(i) ), сложение векторов z = x + y ( z(i) = x(i) + y(i) ), скалярное произведение c = x(T)y ( c = sum(i=1,n)( x(i)y(i) ) ) и покомпоненнтное произведение z = x.*y ( z(i) = x(i)y(i) ). Пятая операция - saxpy - настолько важна для матричных вычислений, что мы дополняем ею наш список основных операций над векторами, несмотря на то, что это делает наш список избыточным. Операция saxpy определяется следующим образом: z = ax + y => z(i) = ax(i) + y(i) Название "saxpy" используется в пакете Linpack, в котором реализованы многие из приведенных в этой книге алгоритмов. Можно воспринимать "saxpy" как мнемонику для выражения "скаляр альфа, умноженный на x, плюс y" (scalar alpha x plus y). 1.1.5. Вычисление скалярных произведений и операций saxpy Для записи алгоритмов мы решили использовать стилизованную версию языка Matlab. Matlab - это элегантная интерактивныя система, идеально приспособленная для вычислений с матрицами. См. Coleman and Van Loan (1988, гл. 4). Наш способ записи вводится постепенно на протяжении этой главы, здесь мы начинаем с функции, вычисляющей скалярное произведение. Алгоритм 1.1.1 (Скалярное произведение.) По n-векторам x и y этот алгоритм вычисляет их скалярное произведение c = x(T)y. function: с = dot(x,y) c = 0 n = length(x) for i = 1:n c = c + x(i)y(i) end end dot В этой процедуре размерность векторов дается фукцией length. Оператор for показывает, что переменная-счетчик i принимает занчения 1, 2, ...., n. Значение скалярного произведения возвращается в с. Скалярное произведение n-векторов требует n умножений и n сложений. Менее формально - это О(n)-операция, в том смысле что объем работы линейно зависит от размерности. Вычисление saxpy - это также O(n)-операция, но она возвращает не скаляр, а вектор: Алгоритм 1.1.2 (Saxpy). По n-векторам x и y и скаляру a этот алгритм вычисляет z = ax + y function: z = saxpy(a,x,y) n = length(x) for i = 1:n z() = ax() + y() end end saxpy 1.1.6. О формализме и стиле изложения Неоходимо подчеркнцть, что приводимые в этой книге алгоритмы, такие как dot и saxpy, лишь заключают в себе важнейшие вычислительные идеи, но никак не являются техническими инструкциями для написания программ. Более того, там, где того требует изложение, записи алгоритмов могут становиться весьма неформальным. Для описания вычислений в этом и в нескольких последующих разделах оказываются удобными функции системы Matlab, но дальше по ходу книги алгоритмы становятся слишком сложными для детализации на "ij"-уровне, и тогда мы прибегаем к менее формальной записи. Для читателя, который по первым главам приобретет интуитивное понимание матричных вычислений, это не должно представлять трудностей. 1.1.7. Умножение матрицы на вектор Пусть A -> R(m x n), и мы хотим вычислить произведение z = Ax, где y -> R(n). Стандартный способ вычисления заключается в последовательном подсчете скалярных произведение z(i) = sum(j=1,n)( a(ij) y(j) ) Это приводит к следующему алгоритму. Алгоритм 1.1.3 (Умножение матрицы на вектор: Версия с доступом по строкам). По A -> R(m x n) и x -> R(n) следующий алгоритм вычисляет z = Ax. function: z = matvec.ij(A,x) m = rows(); n = cols(A) z(1:m) = 0 for i = 1:m for j = 1:n z(i) = z(i) + A(i,j)x(j) end end end matvec.ij Функция rows и cols принимают в качестве аргумента матрицу и возвращают количество ее строк и столбцов соотвественно. Оператор z(1:m) = 0 инициализируюет z как нулевой m x 1-вектор. Альтернативный алгоритм получится, если рассмотреть z = Ax как линейную комбинацию столбцов матрицы A, например [1 2] [1*7 + 2*8] [1] [2] [23] [3 4][7] = [3*7 + 2*8] = 7[3] + 8[4] =[53] [5 6][8] [5*7 + 2*8] [5] [6] [83] Алгоритм 1.1.4 (Умножение матрицы на вектор: Версия с доступом по столбцам). По A -> R(m x n) и R(n) следующий алгоритм вычисляет z = Ax. function: z = matvec.ji (A,x) m = rows(A); n = cols(A); z(1:m) = 0 for j = 1:n for i = 1:m z(i) = z(i) + x(j)A(i,j) end end end matvec.ji Заметим, что цикл по i выполняет операцию saxpy. Мы получили столбцовую версию, переосмыслив на уровне векторов алгоритм матрично-векторного умножения. Другим способом мы могли бы вывести столбцовую версию, поменяв порядок следования циклов в строчном алгоритме. В матричных вычилениях важно осознавать алгебраический смысл модификаций программ, таких как перестановка циклов. 1.1.8. Разбиение матрицы на строки и столбцы Анализ внутренних циклов показывает, что в алгоритме 1.1.3. доступ к элементам матрицы A осуществляется по строкам, а в алгоритме 1.1.4. - по столбцам. Чтобы яснее охарактеризовать эти методы доступа, нам понадобится язык блочных разбиений матриц. С точки зрения строк, матрица представляет собой набор веторов-строк: [a(T)(1)] A -> R(m x n) <=> A = [.......], a(k) -> R(n) [a(T)(m)] Это разбиение матрицы A на строки. Так, разбиение на строки матрицы [1 2] [3 4] [5 6] означает, что мы предпочитаем рассматривать A как набор строк, где a(T)(1) =[1 2] ,a(T)(2) =[3 4],a(T)(3) = [5 6] С учетом разбиения на строки (1.1.1) мы видим, что функция matvec.ij устроена, в сущности, следующим образом: for i = 1:m z(i) = a(T)(i)x end Иначе матрицу можно рассматривать как набор векторов-столбцов: A -> R(m,n) <=> A = [ a(1),....,a(n) ], a(k) -> R(m). Это разбиение матрицы A на стобцы. Для вышеприведенного примера мы положили бы a(1) и a(2) равными первому и второму столбцу матрицы A соотвественно: [1] [2] a(1) = [3] , a(2) =[4] . [5] [6] С учетом разбиения (1.1.2) мы видим, что функция matvec.ji - это использующие saxpy процедура, в которой доступ к элементам матрицы A идет по столбцам: z(1:m) = 0 for j = 1:n z = z + x(j)a(j) end Отметим, что в векторе z накапливается сумма, значение которой обновляется операциями saxpy. 1.1.9. Использование двоеточия Удобным способом сослаться на какой-либо столбец или строку матрицы является использование двоеточия. Обозначим через A(k,:) k-ю строку матрицы A( A -> R(m,n) ): A(k,:) = [a(k1), ....., a(kn)]. Аналогично, обозначим через A(:,k) k-й столбец матрицы A: [a(1k)] A(:,k) = [.....] [a(mk)] Используя эти соглашения, мы можем переписать matvec.ij в виде for i = 1:m for i = 1:m z() = A(i,:)x или z() = dot(A(i,:)x) end end Используя saxpy версия matvec.ji запишется в виде z(1:m) = 0 z(1:m) = 0 for j = 1:n for j = 1:n z = z + x(j)A(:,j)x или z = saxpy( x(j), A(:,j), z ) end end При описании матричного алгоритма важно верно выбрать уровень и стиль описания. Использование двоеточия позволяет исключить из рассмотрения внутренние циклы, тем самым фокусируя внимание на операциях более высокого уровня. Используя функции, мы можем выявить ключевые вычислительные ядра, которые могут отражать хорошие программыне разработки, например dot и saxpy. 1.1.10. Модификация внешним произведением Как мы видели, метод доступа является важным атрибутом матричного алгоритма. По причинам, о которых мы скажем позже, предпочтительными оказываются алгоритмы, выбирающие элементы массивов по столбцам. Так, версия matvec.ji в большинстве случаев предпочитается версии matvec.ij. Вспомним, что единственная разница между этими версиями - порядок следования двух циклов. Важно осознать эту взаимосвязь между порядком следования циклов и методом доступа к данным. Для дальнейшей иллюстрации этого рассмотрим модификацию матрицы A внешним произведением: A <- A + xy(T), A -> R(m,n), x -> R(m), y -> R(n) Здесь <- означает присваивание. Вспоиним, что xy(T) - это просто частный случай произведенияя матриц, так что элементы матрицы A перевычисляются по формуле a(ij) = a(ij) + x(i)y(j), i = 1:m, j = 1:n. Версия "ij" этого вычисления говорит, что надо к каждой строке матрицы A привавить вектор, кратный y(T): for i = 1:m A(i,:) = A() + x(i)y(T) end С другой стороны, версия "ji" for j = 1:n A(:,j) = A(:,j) + y(j)x end выбирает элементы матрицы A по столбцам. Заметим, что обе прецедуры используют saxpy. 1.1.11. Вычисление операции gaxpy Модификация матрицы внешним произведением векторов занимает большое место в традиционных формулировках многих важных матричных алгоритмов. Оказывается, большинство этих алгоритмов можно переформулировать так, что доминирующей становится операция gaxpy. Операция gaxpy - это просто вычисление вида z = y + Ax, x ->R(n), y ->R(m) , A -> R(m,n) Как и "saxpy" , термин "gaxpy" обязан своим происхождением программному обеспечению. Можно рассматривать его как мнемонику для фразы "произвольная матрица A, умноженная на x плюc y" (general Ax plus y). Как будет объяснено в 1.4., формулировки посредством gaxpy (должным образом реализованные), как правило, предпочтительнее формулировок посредством модификации внешним произведением. Это очень важное с вычислительной точки зрения положение, так что gaxpy достойна формального алгоритмического представления: Алгоритм 1.1.5.(Gaxpy). По x -> R(n),y -> R(m) и A -> R(m,n) этот алгоритм вычисляет z = y + Ax function: z = gaxpy(A,x,y) n = cols(A); z = y for j = 1:n z = z + x(j)A(;,j) end end gaxpy Мы организовали вычисления так, что в векторе z накапливается сумма, значение которой обновляется последовательностью операций saxpy. При такой формулировке мы видим, что gaxpy является обобщением saxpy. 1.1.12. Понятие уровня Скалярное произведение и saxpy - это примеры операций "уровня 1". Для операций уровня 1 количество данных и количество арифметических действий линейно завися от размерности операции. Модификация m x n-матрицы внешним произведением или операция gaxpy производять квадратичный объем работы О(mn) над квадратичным количеством данных О(mn). Это примеры операций "уровня 2". Разработка матричных алгоритмов, содержащих много операций высокого уровня, является областью интенсивных исслидований, и мы постоянно будем возвращаться к этой теме в нашей книге. К примеру, для высокой производительности алгоритма решения систем линейных уравнений может потребоваться, чтобы гауссово исключение было организовано через операции gaxpy уровня 2. Для этого понадобится некоторое переосмысление обычного алгоритма, так как в стандартной формулировке гауссово исключение описывается на уровне 1, например, "умножить строку 1 на константу и прибавить результат к строке 2". Кроме того, некоторые алгоритмы можно организовать таким образом, что в них будет много умножений матрицы на матрицу - это операции уровня 3. Операции уровня 3 производят кубический объем работы над квадратичным количеством данных. Если A, B и C - матрицы, то вычисления C = C + AB и C = C + AB - примеры операций уровня 3. Эти состоящие из большого количества действий вычисления могут быть организованы многочисленными спосообами, как мы сейчас прдемонстрируем на примере обычного умножения матриц C = AB. 1.1.13. Умножение матриц Рассмотрим задачу вычисления произведения 2 x 2-матриц C = AB. В терминах скалярных произведений каждый элемент матрицы C вычисляется как скалярное произведение: [1 2][5 6] [1*5 + 2*7 1*6 + 2*8] = [3 4][7 8] [3*5 + 4*7 3*6 + 4*8] В версии, использующей saxpy, каждый столбец матрицы C рассматривается как линейная комбинация столбцов матрицы A: [1 2][5 6] [ [1] [2] [1] [2]] = [5 [ ] + 7[ ]6[ ] + 8[ ]] [3 4][7 8] [ [3] [3] [3] [4]] Наконец, в версии, использующей внешние произведения, C получается как сумма внешних произведений: [1 2][5 6] [1] [2] [ ][ ] =[ ][56] + [ ][78]. [3 4][7 8] [3] [4] Оказывается, что, несмотря на математическую эквивалентность всех трех версий, их реальная производительность может существенно различаться из-за разных доступов к памяти. Мы продолжим обсуждение этого вопроса в 1.4., а пока рассмотрим более детально обозначенные выше три подхода к умножению матриц. Это позволит нам поупражняться в конструировании различных вариантов алгоритмов и записи их с использованием введенных обозначений. 1.1.14. Умножение матриц с использованием скалярных произведений Пусть A -> R(m,r),B -> R(r,n) и мы хотим вычислить C = AB. Стандартная процедура заключается в последовательном вычислениее элементов в порядке слева направо сверху вниз: Алгоритм 1.1.6 (Умножение матриц: версия, использующая скалярное произведение). По A -> R(m,r) и B -> R(r,n) алгоритм вычисляет C = AB. function: C = matmat.ijk(A,B) m = row(A); r = cols(); n = cols(B) C(1:m, 1:n) = 0 for i = 1:m for j = 1:n for k = 1:r C(i,j) = C(i,j) + A(i,k)B(k,j) end end end end matmat.ijk В этой процедуре c(i,j) вычисляется как скалярное произведение i-й строки матрицы и j-го столбца матрицы B. На языуе блочных матриц, если мы имеем разбиения [a(1)(T)] A =[. ], a(k) -> R(r) [a(m)(T)] и B = [b(1),.......b(n)], b(k) ->R(r) то алгоритм 1.1.6 описывается на уровне 1 как for i = 1:m for j = 1:n c(ij) = a(i)(T)b(j) end end Заметим, что задача цикла по j- вычислить i-ю строку матрицы C. Мы можем подчеркнуть это, записав for i = 1:m c(i)(T) = a(i)(T)B end где c(i)(T)- i-я строка C. Это описание алгоритма в терминах матрично-векторных произведений можно рассматривать как описание уровня 2. Конечно, на самом верхнем уровне останется запись C = AB, не содержащая ни циклов, ни индексов. Правильный выбор уровня описания важен как для изложения теоретических результатов, так и для разработки программного обеспечения. 1.1.15. Умножение матриц с использованием gaxpy Предположим, что A, B и C разбиты на столбцы: 25 A = [a(1),....a(r)], a(j) -> R(m)? B = [b(1),....b(n)], b(j) -> R(r)? C = [c(1),....c(n)], c(j) -> R(m)? и мы желаем вычислить C = AB. Поскольку c(j) = sum(k,1,r)( b(k,j)a(k) ), j = 1:n, мы видим, что каждый столбец C является линейной комбинацией столбцов A. Вычисление этих сумм можно организовать как последовательность операций gaxpy: Алгоритм 1.1.7. (Умножение матриц: gaxpy-версия). По A -> R(m,r) B -> R(r,n) и следующий алгоритм вычисляет C = AB. function C = matmat.jki(A,B) m = rows(A); n = cols(B); r. = cols(A) C(1:m, 1:n) = 0 for j = 1:n for k = 1:r for i = 1:m C(i,j) = C(i,j) + A(i,k)B(k,j) end end end end matmat.jki Мы вывели matmat.jki, разбив по столбцам матрицы в равенстве C = AB, но можно было получить его также изменением порядка циклов в функции matmat.ijk. Ясно, что "jki"-алгоритм состоит из операций gaxpy, поскольку его можно записать в виде C(1:m, 1:n) = 0 for j = 1:n C(:,j) = gaxpy(A, B(:,j),(:,j)) end 1.1.16. Умножение матриц с использованием внешних произведений Давайте проделаем перестановку циклов еше раз, а затем проанализируем ее смысл в терминах линейной алгебры. Алгоритм 1.1.8(Умножение матриц: версия, использующая внешние произведения). По A -> R(m,r) и B -> R(r,n) алогритм вычисляет C = AB. function: C = matmat,kji(A,B) m = rows(A); n = cols(B); r = cols(A) C(1:m, 1:n) = 0 for k = 1:r for j = 1:n for i = 1:m C(i,j) = C(i,j) + A(i,k)B(k,j) end end end end matmat.kji В этой "kji"-формулировке происходит вычисление всех элементов c(ij) одновременно. Для данного k два внутренних цикла осуществляют модификацию матрицы C внешним произвадением C <- C + a(k)b(k)(T), где A = [a(1),....,a(r)], a(k) -> R(m) и [b(1)(T)] и Внутренний цикл в функции matmat.kji выполняет операцию saxpy. Именно, к j-му столбцу матрицы C прибавляется вектор, кратный a(k). 1.1.17. Перестановка циклов Двойной цикл в алгоритме умножения матриы на вектор может быть переупорядочен 2! = 2 способами. Соотвественно имеется 3! = 6 способов упорядочивания тройного цикла в алгоритме умножения матрицы на матрицу. Три варианта уже были в деталях изложены выше: ijk, jki, и kji. Каждый из шести возможных вариантов характеризуется своим способом доступа к данным и доминирующей операции (скалярное произведение, saxpy). Эти сведения собраны в табл. 1.1.1. Какой же вариатн предпочтительнее? Ответ на этот вопрос зависит от архитектуры компьютера, на котором будет реализован алгоритм. Подробнее это обсуждается в 1.4. Таблица 1.1.1. Матричное умножение; упорядочение и свойсва Порядок Внутренний Средний Внутренний цикл цикла цикл цикл Доступ к данным =============================================================== ijk скалярное вектор x матрица A на строку,B на столбец произведение --------------------------------------------------------------- jik скалярное матрица x вектор A на строку,B на стобец --------------------------------------------------------------- ikj saxpy строчное gaxpy B на строку ---------------------------------------------------------------- jki saxpy столбцовое gaxpy B на столбец ---------------------------------------------------------------- kij saxpy строка внешнего про- B на строку изведения ---------------------------------------------------------------- kji saxpy столбец внешнего A на столбец произведения ================================================================ 1.1.18. О матричных соотношениях Изучая умножение матриц с использованием внешних произведений, мы по существу установили соотношение AB = sum(k=1,r)( a(k)b(k)(T) ) где a(k) и b(k) определены в (1.1.3) и (1.1.4). 27 В последующих главах мы получим очень много матричных соотношений. Иногда они выводятся алгоритмически, как это было в приведенном выше примере с внешним произведением, а в других случаях их приходится доказывать на уровне "ij"-компонент. Как пример этого стиля доказательств приведем важный результат, характеризующий транспонирование произведения матриц. Теорма 1.1.1. Если A -> R(m,r) и B -> R(r,n), то (AB)(T) = B(T)A(T). Доказательство. Если C = (AB)(T), то c(ij) = [(AB)(T)](ij) = (AB)(ji) = sum(k=1,r)( a(jk)b(ki)). С другой стороны, если D= B(T)A(T), то d(ij) = [B(T)A(T)](ij) = sum(k=1,r)( B(T)(ik)(A)(T)(kj) ) = sum(k=1,r)( b(ki)a(jk) ), так что C = D. Как свидетельствует это доказательство, погружение на "ijk"-уровень отнюдь не тождественно глубокому прокникновению в материал. Однако иногда это единственный способ получить алгебраический результат высокого уровня. 1.1.19. Об описании алгоритмов На протяжении этой книги нам понаюобится описывать алгоритмы на разных уровнях детализации. В этом разделе все алгоритмы были записаны формально, как функции. В качестве примера более вольного стиля записи алгоритмов приведены процедуру, вычисляющую матричное произведение . Алгоритм 1.1.9. По двум m x n-матрицам A и B алгоритм вычисляет произведение C = A(T)B. for i = 1:n for j = 1:n C(ij) = A(1:m,i)(T)B(1:m,j) 1.1.20. Комплексные матрицы Наше внимание в этой книге сфокусировано на вычислениях с вещественными матрицами - для упрощения обозначений, а также ввиду того, что большинство реальных задач имеет дело с вещественныи данными. Однако в последующих разделах мы будем работать и с комплексными матрицами - там, где это уместно. Для начала введем некоторые обозначения. Векторное пространство комплексных m x n-матриц обозначается через C(m,n). Умножение на скаляр, сложение и перемножение комплексных матриц такое же, как и в вещественном случае. Однако транспонирование превращается в сопряженное транспонирование: C = A(H) -> c(ij) = (a-)(ji) Векторное пространство комплексных n-векторов обозначается C(n). Скалярное произведение комплексных n-векторов вычисляется по правилу s = x(H)y = sum(i=1,n)( (x-)(i)y(i) ) Наконец, если A = B + iC -> C(m,n), то мы будем обозначать вешественную и мнимую части A через Re(A) = B и Im(A) = C соотвественно. Замечания и литература к $ 1.1 Подчеркнем еще раз, что приводимые в этой книге алгоритмы не являются готовым программным продуктом. Если подходить к делу всерьез, то написание качественной программы на основе любого из наших описаний алгоритмов требует дилтельной и напряженной работы. Даже реализация процедур низкого уровня, таких как saxpy и gaxpy, требует времени и внимания. Подробнее эти вопросы рассмотрены в статьях: Dongarra J., DuCroz J., Hammarling S., and Hanson R.J.(1988). "An Ehtended Set of Fortran Basic Linear Algebra Subprograms", ACM Trans. Math. Soft. 14, 1-17. Dongarra J., DuCroz J., Hammarling S., and Hanson R.J.(1988). "Algorithm 656 An Ehtended Set of Fortran Basic Linear Algebra Subprograms: Model Impleventation and Test Prograns", ACM Trans. Math. Soft. 14, 18-32. Dongarra J., DuCroz J., Duff I.S., and Hammarling S.(1988). "A Set of Level 3 Basic Linear Algebra Subprograms", Argonne National Laboratory Report, ANL-MCS-TM-88. Lawson C.L., Hanson R.J., Kincaid D.R., and Krogh F.T.(1979a). "Basic Linear Algebra Subprogram for FORTRAN Usage", ACM Trans. Math. Soft. 5, 308-23. Lawson C.L., Hanson R.J., Kincaid D.R., and Krogh F.T.(1979b). "Algorithm 539, Basic Linear Algebra Subprogram for FORTRAN Usage", ACM Trans. Math. Soft. 5, 324-25. Мы также рекомендуем книгу Rice J.R.(1981). Matrix Computations and Mathematical Software, Academib Press, New York. [Имеется русский перевод: Райс Дж. Матричные вычисления и математическое обеспечение.-М.: Мир, 1984.] Влияние различных упорядочений цеклов на производительность детально рассмотрено в великолепном обзоре Dongarra J.J., Gustavson F.G., and Karp A.(1984). "Implementing Linear Algebra Algorithms for Dense Matrices on a Vector Pipeline Machine", SIAM Review 26, 91-112. 1.2 Учет структуры матрицы 29 Эффективность заданного матричного алгоритма зависит от многих характеристик. В этом разделе мы рассмотрим наиболее очевидню из них - количество арифметических операций и объем памяти, тебуемые процедурой. Менее очевидны такие характеристики, как способ доступа к памяти, шаг выборки и еще целый сонм издержек, возникающих при пересылке данных. Эти вопросы рассматриываются, в частности, в 1.4. Мы продолжаем обкатывать новые ключевые идеи на примерах умножения матрицы на матрицу и матрицу на вектор. Для деманстрации учета структуры матриц мы избрали свойства ленточнотси и симметричности. Ленточные матрицы содержат много нулевых элементов, поэтому неудивительно, что умножение ленточных матриц позволяет уменьшить по сравнению в общим случаем и арифметические затраты, и объем требуемой памяти. В этом контекстве мы обсуждаем структуры данных и понятие "флопа". Другой пример учета особенностей структуры доставляют симметричные матрицы. Решение систем линейных уравнений и задач на собственные значения с симметричными матрицами играет видную роль в матричных вычислениях, поэтому важно приобрести хороший навык работы с ними. Мы заканчиваем этот раздел замечаниями о записи результатов вычислений на место входных данных - это еще один прием для контроля объема хранимо информации. 1.2.1. Ленточные матрицы и обозначения для них 1.2.2. Действия с диагональными матрицами 1.2.3. Умножение треугольных матриц 30 16 1.2.4. Флопы Очевидно, умножение верхних треугольных матриц требует меньшего количеста арифметических операция, чем умножение произвольных матриц. Один из способов количественной характиристики разницы дает понятие флопа. Флоп - это одна операция с плавающей точкой (floating point operation). Скалярное произведение или операция saxpy размерности n сдержит 2*n флопов, поскольку каждая из них состоит из n умножение и n сложений. 1.2.5. Еще раз об использовании дваеточия Скалярное произведение, вычисляемое циклом по k в алгоритме 1.2.1, можно записать компактно, если обобщить введенное в 1.1.9 использование двоеточия. Пусть A -> R(m,n) и целые числа p, q и r удовлетворяют неравенствам 1 <= p <= q <= n и 1 <= r <= m. Тогда определим A(r,p:q) = [a(rp),.....a(rq)] -> R(1,(q - p +1)). Аналогично, если 1 <= p <= q <= m и 1 <= c <= n, то [a(pc)] A(p:q,c) = [ . ] -> R(q - p + 1) [a(qc)] С использованием этих обозначений алгоритм 1.2.1. перепишется в виде: C(1:n, 1:n) = 0 for i = 1:n for j = 1:n C(i,j) = C(i,j) + A(i,i:j)B(i:j,j) end end Упомянем еще одно добавление к нашим обозначениям. При использованиее двоеточия допускаются отрицательные приращения. Так, если x и y - n-векторы, то s = x(T)y(n:-1:1) есть сумма s = SUM(i=1,n)x(i)y(n - i + 1), т.е. свертка векторов x и y. 1.2.6. Хранение ленточных матриц 1.2.7. Симметрии 1.2.8. Хранение по диагоналям 1.2.9. Запись поверх входных данных 1.3 Блочные матрицы и алгоритмы 36 Умение непринужденно работать с блочными матрицами крайне необходимо для матричных вычислений. Использование блочных обозначений упрощает вывод многих важных алгоритмов. Кроме того, все более возрастает роль блочных алгоритмов в вычислениях на суперкомьютерах. Здесь под блочныхми алгоритмами мы подразумеваем по существу такие алгоритмы, которые содержат много операций умножения матрицы на матрицу. Следует ожидать, что эти алгоритмы будут более эффективны по сравнению с теми, которые работают на скалярном уровне, поскольку на каждый обмен данных будет приходится большее количество вычислений. К примеру, умножение k x k-матрицы занимает 2k**3 флопов, вовлекая лишь 2k**2 элементов данных. В 1.4. мы увидим, что при больших k издержки обмена становятся относительно менее существенными. 1.3.1. Обозначения для блочных матриц 1.3.2. Операции с блочными матрицами 1.3.3. Обозначения для подматриц 1.3.4. Умножение блочной матрицы на вектор 1.3.5. Перемножение блочных матриц 1.3.6. Важный частный случай 1.3.7. Структуры данных для блочных матриц 1.3.8. Умножение матриц, основанное на принципе "разделяй и властвуй" 45 (23) 1.4 Некоторые аспекты векторно-конвейерных вычислений 45 Действия с матрицами в основном слагаются из скалярных произведений и операций. Векторно-конвейерные компьютеры способны очень быстро выполнять подобные операции благодаря специальной аппаратной поддержке, учитывающей то обстоятельство, что векторная операция - это очень регулярно устроенная последовательность скалярных операций. Будет ли при использовании такого компьютера достигнута высокая производительность, это зависит от длины векторных операндов, а также ряда факторов, относящихся к пересылкам данных, таких как шаг выборки, количество векторных загрузок/записей, степень повторного использования вектора. Знакомство с этими вопросами и понимание их роли чрезвычайно полезны. Мы не будем пытаться построисть всеобъемлющую модель векторно-конвейерных вычислений, позволяющую предсказывать производительность алгоритмов. Наша цель - дать читателям почувствовать основные моменты, которыми следует руководствоваться при создании эффективных программ для векторно-конвейерных компьютеров. Мы не ориентируемся на какую-то конкретную машину; заинтересованных читателей мы отсылаем к литерауре, которая содержит очень много конкретных исследований такого рода. 1.4.1. Конвейеризация арифметических операций 1.4.2. Векторные операции 1.4.3. Длина вектора 1.4.4. Оптимизация кода с учетом длины векторов 1.4.5. Многоуровневая память 1.4.6. Единичный шаг выборки 1.4.7. Роль структур данных 1.4.8. Операции , внешние произведения и вектоные обмены 1.4.9. Блочные алгоритмы, кэш-память и повторное использование данных 1.4.10. Резюме Глава 2. Матричный анализ 56 Вывод алгоритмов для матричных вычислений и их анализ требуют хорошего знакомства с некоторыми вопросами линейной алгебры. Обзор ряда основных понятий линейной алгебры дается в . В и рассмотрены векторные и матричные нормы. В мы строим модель арифметики с конечной точностью и знакомим читателя с методами матричного анализа, потребными для количественной оценки влияния ошибок округления. Следующие два раздела посвящены ортогональности, которая играет в матричных вычислениях заметную роль. Две разновидности ортогональных приведений - сингулярное разложение и CS-разложение - дадут нам необходимую теоретическую основу для введения важных понятий ранга и расстояния между подпространствами. Эти вопросы излагаются в и . Наконец, в последнем разделе мы изучаем поведение решения линейной системы при возмущении матрицы и вектора. При этом вводится важное понятие числа обусловленности матрицы. 2.1 Основные сведения из линейной алгебры 56 Этот раздел представляет собой беглый обзор линейной алгебры. Читатели которые пожелают ознакомится с более подробным изложением предмета, мы рекомендуем обратиться к ссылками в конце этого раздела. 2.1.1. Линейная независимость, подпространства, базис и размерность 2.1.2. Область значений, ядро и ранг матрицы 2.1.3. Обратная матрицы 2.1.4. Детерминант 2.1.5. Дифференцироание Замечания и литература к 2.1. Среди множенства вводных курсов линейной алгебры только некоторые обеспечивают ничинающего изучать матричные вычисления необходимыми ему материалом. Особенно полезными мы считаем следующие книги: Halmos P.R.(1958). Finite Dimensional Vector Spaces, 2nd ed., Van Nostrand-Reinhold, Princeton. [Имеется русский перевод: Халмош П.Р. Конечномерные векторные пространства.-М.: Физматгиз, 1960.] Leon S.J.(1980). Linear Algebra with Applications. Macmillan, New York. Noble B. and Daniel J.W.(1977). Applied Linear Algebra, Precntice-Hall, Englewood Cliffs. Strang G.(1988). Linear Algebra and Its Applications, 2rd edition, Harcourt, Brace, Jovanovich, San Diego. Более энциклопедическое изложение можно найти в Bellman R.(1970). Introduction to Matrix Analysis, 2nd ed., McGraw-Hill, New York. [Имеется русский перевод первого издания: Беллман Р. Введение в теорию матриц.-М.:Наука, 1969.] Gantmacher F.R.(1959). The Theory of Matrices, vols. 1 and 3, Chelsea, New York. [Имеется русский оригинал: Гантмахер Ф.Р. Теория матриц. 4 изд.-М.:Наука, 1989.] Householder A.S.(1964). The Theory of Matrices in Numerical Analysis, Ginn (Blaisdell), Boston. Стюарт (IMC, гл.1) также дает великолепный обзовр матричной алгебры. 2.2 Векторыне нормы 59 Нормы в векторных пространствах служат той же цели, что и взятие модуля на вещественной прямой: они позволяют измерять расстояние. Точнее, R(**n) вместе с введенной на R(**n) нормой определяет матрическое пространство. Поэтому при работе с векторами и векторзначными функциями мы будем располагать привычными понятоиями окрестности, открытого множетсва, сходимости и непрерывности. 2.2.1. Определения 2.2.2. Некоторые свойства векторных норм 2.2.3. Абсолютная и относительная погрешности 2.2.4. Сходимость 2.3 Матричные нормы 61 Матричные нормы часто требуются при анализе матричных алгоритмов. Например, программа решения систем линейных уравнений может давать некачественный результат, если матрица коэффициентов "почти вырожденная". Для количественной характеристики близости к вырожденности нам нужно умет измерять расстояния в пространстве матриц. Такую возможность дают матричные нормы. 2.3.1. Определения 2.3.2. Некоторые свойства матричных норм 2.3.3. Матричная 2-норма 2.3.4. Возмущения и обратная матрица 2.4 Матричные вычилсления с конечной точностью 65 Ошибки округления - это одна из тех вещей, которые делают матричные вычисления нетривиальной и интересной областью. В этом разделе мы построим модель арифметики с плавающей точкой и применим ее для получения оценок ошибок при вычислении с плавающей точкой скалярных произведений, операций saxpy, произведений матрицы на вектор и матрицы на матрицу. 2.4.1. Числа с плавающей точкой 2.4.2. Модель арифметики с плавающей точкой 2.4.3. Потеря точности 2.4.4. Обозначение абсолютный величины 2.4.5. Ошибки округления в скалярных произведениях 2.4.6. Альтернативные методы количественной оценки ошибок 2.4.7. Вычисление скалярных произведений с накоплением 2.4.8. Ошибки округления в других основных матричных вычислениях 2.4.9. Прямой и обратный анализ ошибок 2.4.10. Ошибки округления в алгоритме Штрассена Замечания и литература к 2.4 Наиболее полное изложение анализа ошибок округления содержится у Уилкинсона (AEP, гл. 3). Первосходное изложение имеется также у Форсайта и Молера (SLE, pp. 87-97) и Стьюарта (IMC, pp. 69-82). Для общего знакомства с ролью ошибок округления мы рекомендуем гл. 2 книги Forsythe G.E., Malcolm M.A., and Moler C.B.(1977). Computer Methods for Mathematical Computations, [Имеется русский перевод: Форсайт Дж., Малькольм М., Молер К. Машинные методы математических вычислений.-М.:Мир, 1980] и ее переработанного издания Kahaner D., Moler C.B., and Nash S. (1988). Numerical Methods and Software. Мы уверенно раекомендуем классическое руководство Wilkinson J.H.(1963). Rounding Errors in Algebraic Processes. Философско-истрорический обзор анализа ошибок округления содержистя в фон-неймановской лекции Уилкинсона, опубликованной в Wilkinsom J.H.(1971). "Modern Error Analysis", SIAM Review 13, 548-568. В этой статье дано критическое рассмотрение ранних работ по анализу ошибок, принадлежащих фон Нейману и Голдстайну, Тьюрингу и Гивенсу. Более современные течения в анализе ошибок включают интервальный анализ, построение статистических моделей ошибок округления и автоматизацию самого процесса анализа. Смотрите статьи Hull T.E. and Swensen J.R. (1966). "Tests of Probabilistic Models for Propagation of Roundoff Errors", Comm. ACM 9, 108-113. Larson J. and Sameth A.(1978). "Efficient Calculation of the Effects of Roundoff Errors", ACM Trans. Math. Soft. 4, 228-236. Miller W. and Spooner D.(1978). "Software for Roundoff Analysis, II", ACM Trans, Math. Soft. 4, 369-390. Yohe J.M.(1979). "Sofrware for Interval Arithmetic: A Reasonable Portable Package", ACM Trans. Math. Soft. 5, 50-63. Всякому, кто всерьез занимается разработкой программного обеспечения, необходимо глубокое понимание специфики вычислений с плавающей точкой. Неплохим отправным пунктом для приобретения знаний в этой области станет ознакомление со стандартом IEEE арифметики с плавающей точкой по работе Stevenson D.(1981). "A Proposed Standard for Binary Floating Point Arithmetic", Computer 14 (March), 51-62. Желательные свойства системы чисел с плавающей точкой по-прежнему интенсивно разрабатываются. Смотрите статьи Demmel J.W.(1984). "А Proposed Standard for Binary Floating Point Arithmetic", Computer 14 (March), 51-62. Kulish U.W. and Mirenker W.L.(1986). "The Arithmetic of the Digital Computer", SIAM Review 28, 1-40. Даже для простейших задач разработка высококачественного программного обеспечения сопряжена с неимоверным количеством тонкостей и нюансов. Хороший пример - разработка подпрограммы для вычисления 2-норм: Blue J.M.(1978). "A Portable FORTRAN program to Find the Euclidian Norm fo a Vector", ACM Trans. Math. Soft. 4, 15-23. Анализ метода Штрассена и других "быстрых" алгоритмов линейной алегебры можно найти в Brent R.P.(1970). "Error Analysis of Algorithms for Matrix Multiplication and Triangular Decomposition Using Winograd's Identity", Numer. Math. 16, 145 - 156. Miller W.(1975). "Computational Complexity and Numerical Stability", SIAM J. Computing 4, 97-107. 2.5 Ортогональность и сингулярное разложение 73 Ортогональность играет в вычислениях с матрицами выдающуюся роль. Введя несколько определений, мы докажем очень полезную теорему о сингулярном разложении (SVD). Помимо всего прочего, SVD позволяет по-умному подойти к проблеме ранга матрицы. Понятие ранга, совершенно ясное в условиях точных вычислений, становится скользким в присутствии ошибок округления и прогрешностей в исходных данных. Сингулярное разложение позволяет ввести практическое понятие численного ранга. 2.5.1. Ортогональность 2.5.2. Нормы и ортогональные преобразования 2.5.3. Сингулярное разложение 74 2.5.4. Неполнота ранга и SVD 2.6 Проекции и CS-разложение 77 Если в результате вычисления должна получиться матрица или вектор, то для оценки точности ответа или для измерения того, насколько успешно продвигается итерационный процесс, удобно применять нормы. Если же необходимо вычислить целое подпространство, то для получения аналогичныхх оценок нам нужно умет давать количественное выражение рассторяния между подпространствами. Очень важную роль при этом играют ортогональные проекторы. Рассмотрев элементарные понятия, мы переходим к обсуждению CS-разложения. Это похожее на SVD разложение оказывается очень кстати, когда нужно сравнить пару подпространтсв. Мы начинаем с понятия ортогонального проектора. 2.6.1. Ортогональные проекторы 2.6.2. Проекторы, связанные с SVD 2.6.3. Расстояния, связанные в подпространствами 2.7 Чувствительность квадратных систем к возмущениям 80 В этом разделе мы применяем разработынный в предыдущих параграфах инструментарий к анализу линейной системы Ax = b, где A -> R(n,n) - невырожденная матрица и b -> R(n). Наша цель - исследовать, как возмущения в A и b влияют на решение x. 2.7.1. SVD-анализ 2.7.2. Обусловленность 2.7.3. Определители и близость к вырожденности 2.7.4. Точная оценка по норме 2.7.5. Несколько точных покомпонентных оценок Глава 3. Линейные системы общего вида 87 Проблема решения линейной системы Ax = b является центральной в научных вычислениях. В этой главе мы остановимся на методике исключения Гаусса - методе, который используют, когда матрица A квадратная, плотная и без специфики. Если A не удовлетворяет этим условиям, то представляют интерес алгоритмы из гл. 4, 5 и 10. Параллельные методы решения системы Ax = b обсуждаются в гл. 6. Мы приходим к методу исключения Гаусса, обсуждая в 3.1 ту легкость, с которой можно решать треугольные системы. Приведение системы общего вида к треугольной форме при помощи преобразований Гаусса описывается ниже в 3.2., где выводится язык матричных разложений. К несчастью, полученный метод ведет себя очень плохо на нетривиальном классе задач. Наш анализ ошибок округления в 3.3. выявляет трудности и подготавливает 3.4., где вводится концепция ведущих елементов. В последнем разделе мы предлагаем некоторые замечания по посоду важных для практики вопросов, свзязанных с масштабированием, итерационным уточнением и оценкой обусловленности. 3.1 Треугольные системы 87 Традиционные методы разложений для линейных систем включает в себя приведение исходной квадратной системы к треугольной системе, которая имеет такое же решение. Данный раздел посвящен решению треугольных систем 3.1.1. Прямая подстановка 3.1.2. Обратная подстановка 3.1.3. Столбцовые версии 3.1.4. Случай нескольких правых частей 3.1.5. Доля флопов 3 уровня 3.1.6. Решение неквадратных треугольных систем 3.1.7. Унитреугольные системы 3.1.8. Алгебраические свойства треугольных матриц 3.2 LU-разложение 92 Как мы только что видели, треугольные системы решаются "легко". Идея исключения Гаусса - это преобразование данной системы Ax = b в эквивалентную треугольную систему. Преобразование достигается составлением соответсвующих линейных комбинаций уравнений. Например, в системе 3x + 5y = 9, 6x + 7y = 4, умножая первую строку на 2 и вычитая ее из второй, мы получим 3x + 5y = 9, -3y = -14, Это и есть исключение Гаусса при n = 2. Наша цель в данном разделе - дать полное описание этой важной процедуры, причем описать ее выполнение на языке матричных разложений. Данный пример показывает, что алгоритм вычисляет нижнюю унитреугольную матрицу L и верхнюю треугольную матрицу U так, что A = LU, т.е. [3 5] [1 0][3 5] = [6 7] [2 1][0 -3] Решение для исходной задачи Ax = b находится посредством последовательного решения двух треугольных систем: Ly = b, Ux = u -> Ax = LUx = Ly = b. LU-разложение - это "высокий уровень" алгебраического описания исключения Гаусса. Представление результата матричного алгоритма на "языке" матричных разложений полезно. Оно облегчает обобщение и проясняет связь между алгоритмами, которые могут казаться очень разными на скалярном уровне. 3.2.1. Матрица преобразования Гаусса 3.2.2. Применение матриц преобразования Гаусса 3.2.3. Свойства ошибок округления в преобразованиях Гаусса 3.2.4. Приседение к верхнему треугольному виду 3.2.5. LU-разложение 3.2.6. Несколько практических замечаний 3.2.7. Где хранить матрицк L? 3.2.8. Решение линейной системы 3.2.9. Gaxpy-верси LU-разложения 3.2.10. Модификация Краута-Дулитла 3.2.11. Блочное LU-разложение 3.2.12. LU-разложение прямоугольной матрицы 3.2.13. Несостоятельность метода 3.3 Анализ ошибок округления в методе исключения Гаусса 102 В первых двух разделах мы оценим влияние ошибок округления, когда алгориты используются для решения линейной системы Ax = b. Прежде чем приступить к этому анализу, полезно рассмотреть ситуацию, близкую к идеальной, при которой не возникает ошибок округления в процессе решения, а ошибки возникают при задании A и b. Таким образом, если fl(b) = b + e, а хранимая матрица fl(A) = A + E невырожддена, то мы допускаем, что вычисленное решение x(^) удовлетворяет соотношению (A + E)x(^) = (b + l), ||E|| <= u||A||, ||e|| <= u||b||, (3.3.1) т.е. x(^) является решением системы, близкой к точной. Более того, если uk(A) <= 1/2 (например), то, используя теорему 2.7.1, можно показать, что ||x - x(^)||/||x|| <= 4uk(A). (3.3.2) Оценки (3.3.1) и (3.3.2) являются "наилучшими" оценками по норме. В общем случае не существует анализа ошибок округления в inf-норме для решения линейной системы, который требовал бы хранения A и b и позволял получить более тонкую оценку. Как следствие этого, мы не врпаве критиковать алгоритм за нахождение неточного решения x(^), если матрица плохо обусловлена относительно машинной точности, например uk(A) = 1. 3.3.1. Ошибки в LU-разложении 3.3.2. Решение треугольных систем с приближенными треугольными матрицами 3.4 Выбор ведущего элемента 106 Приведенный в предыдущем разделе анализ показывает, что мы должны быть уверены в отсутствии больших элементов в вычисленных треугольных множителях L и U. Пример [0.0001 1] [ 1 0][0.0001 1 ] A = [ 1 1] = [10000 1][ 0 -9999] = LU показывает источник возникших неприятностей: относительно малый ведущий элемент. Перестановкой строк от этих сложностей можно избавиться. В нашем примере, если P - это матрица перестановок вида [0 1] P = [1 0], то [ 1 1] [ 1 0][1 1 ] PA = [0.0001 1] = [0.0001 1][0 0.9999] = LU Сейчас треугольные множители содержат достаточно малые элементы. В этом разделе мы покажем, как определить перестановочный вариант матрицы A, чтобы добиться достаточно устойчивого LU-разложения. Существует несколько способов достижения цели, и все оин соответсвуют разным стратегиям выбора ведущего элемента. Мы остановимя на стретегиях частичного и полного выбора ведущего элемента. Мы обсудим эффективные реализации этих стратегий и их свойства. А начнем с обсуждения применения матрицы перестановок. 3.4.1. Перестановочные матрицы 3.4.2. Частичный вабор ведущего элемента: общая идея 3.4.3. Детали стратегии частичного выбора 3.4.4. Где хранить матрицу L? 3.4.5. Gaxpy-версия 3.4.6. Анализ ошибок округления 3.4.7. Метод блочного исключения Гаусса 3.4.8. Полный выбор ведущего элемента 3.4.9. Замечания по стратегии с полным выбором ведущего элемента 3.4.10. Отказ от выбора ведущего элемента 3.4.11. Некоторые приложения 3.5 Уточнение и оценивание точности 119 Пусть для решения n x n системы Ax = b используется метод исключения Гаусса с частичным выбором ведущего элемента. Предположим, применяется t-разрядная плавающая арифметика с основанием beta. Уравнение (3.4.4) гарантирует, что если фактор роста умеренный, то вычисленное решение x(^) удовлетворяет соотношению (A + E)x(^) = b, ||E|| = u||A||, u = 1/beta**t В этом разделе мы исследуем практические разветвления этого результата. Сначала сделаем акцент на различии, которое следует делать между величиной невязки и точностью. Потом продолжим обсуждение масштабирования, итерационного уточнения и оценки обусловленности. Прежде чем начать сделаем два замечания по поводу обозначений. Всюду используется бесконечная норма, так как это удобно для анализа ошибок округления и для практической оценки ошибки. Кроме того, где бы в этом разделе мы ни употребили выражение "исключение Гаусса", мы всегда будем имет в виду исключение Гаусса с некоторой устойчивой стратегией выбора ведущего элемента, такой, как частичный выбор. 3.5.1. Большая невяжка дает плохую точность 3.5.2. Масштабирование 3.5.3. Итерационное уточнение 3.5.4. Оценка обусловленности Глава 4. Линейные системы специального вида 127 Основной принцип численного анализа состоит в том, что при решении задач всякий раз надо использовать их специфику. Можно надеятся, что в численных методах ленейной алгебры алгоритмы решения задач с матрицами общего вида приспособлены к использованию таких свойств как симметричность, определенность, разреженность. Это и будет центральной темой данной главы, в которой наша главная цель - предложить специальные алгоритмы для вычисления специальных вариантов LU-разложения. Мы начнем с установления зависимости между множителями L и U, когда матрица A является симметричной. Это достигается исследованием LDMt-разложения в 4.1. Потом мы займемся важным случаем, когда матрица A является одновременно симметричной и положительно определенной, и получим в 4.2. устойчивое разложение Холецкого. В этом разделе мы также исследуем несимметричные положительно определенные системы. В 4.3. обсуждается ленточный вариант исключения Гаусса и другие методы разложения. Потом будет исследована интересная ситуация, когда A симметричная, но неопределенная. Наша трактовка этой задачи в 4.4. выявляет двойственность стоящей перед исследователем проблемы выбора ведущего элемента. Мы предпочитаем выбирать ведущий элемент так, чтобы добиться устойчивости, и не обращаем внимания, что при этом может быть нарушена специфика. К счастью, в задаче с симметричной неопределенной матрицей имеется удачный способ разрешения этого конфликта. Любые блочные ленточные матрицы являются также просто ленточными, и поэтому к ним могут применяться методы из 4.3.. Однако бывают ситуации, когда такой подход неприемлем. Для иллюстрации мы рассмотрим в 4.5. случай блочной трехдиагональной системы. В последующих параграфах мы исследуем несколько очень интересных O(nn) алгоритмов, которые могут быть использованы для решения систем Вандермонда и тёплицевых систем. 4.1 Разложения типа LDM(T) и LDL(T) 127 Мы хотим создать метод, использующий специфику при решении задачи Ax = b. Для этого установим возможность LU-разложения матрицы A в виде произведения трех матриц LDM(T), где D - диагональная, а L и M - нижние треугольные. Когда такое разложение получено, решение Ax = b может быть найдено за O(n**2) флопов посредством решения систем Ly = b (прямая подстановка), Dz = y и M(T)x = z (обратная подстановка). Исследование LDM(T)-разложения позволяет для симметричного случая установить следующий факт: если A = A(T), то L = M и работа, связанная с разложением, составляет половину от того, что требуется для исключения Гаусса. Проблема выбора ведущего элемента поднимается в последующих разделах. 4.1.1. LDMt-разложение 4.1.2. Симметрия и LDLt-разложение 4.2 Положительно определенные системы 132 Матрица A -> R(n,n) является положительено определенной, если x(T)Ax > 0 для всех ненулевых векторов x -> R(n). Положительно определенные системы составляю один из наиболее важных классов задач Ax = b со спецификой. Рассмотрим 2 x 2-симметричный случай. Если [a(11) a(12)] A = [a(21) a(22)] является положительно определенной, то x = (1,0)(T) -> x(T)Ax = a(11) > 0, x = (0,1)(T) -> x(T)Ax = a(22) > 0, x = (1,1)(T) -> x(T)Ax = a(11) + 2a(12) + a(22) > 0, x = (1,-1)(T) -> x(T)Ax = a(11) - 2a(12) + a(22) > 0, Последние два уравнения приводят к неравенству |a(12)| <= (a() + a(22))/2. Из этих результатов мы видми, что наибольший элемента матрицы A лежит на диагонали и он положителен. Это веро и в общем случае. положительно определенная матрица имеет весомую даигональ. Масса диагонили не столь заметна, как в случае с диагональным доиинированием (ср. 3.4.9), но она дает аналогичный эффект - устраняет необходимость выбора ведущего элемента. Мы начнем с нескольких замечаний о свойстве полоижтельной определеннолсти и о том, что она влечен за собой в несимметричном случае при выборе ведущего элемента. Потом мы остановимся на эффективной реализации процедуры Холецкого, которая может быть использована для надежного разложения симметричной положительно определенной матрицы A. Исследуются gaxpy-версия, версия с внешним произведением и блочная версия. Раздел заканчивается несколькими замечаниями о неотрицательно определенном случае. 4.2.1. Положительная определенность 4.2.2. Несимметричные положительно определенные системы 4.2.3. Симметричные полжительно определенные системы 4.2.4. Gaxpy-версия разложения Холецкого 4.2.5. Метод Холецкмого с внешним произведением 4.2.6. Блочный метод Холецкого со скалярным произведением 4.2.7. Устрйчивость процесса Холецкого 4.2.8. Неотрицательно определенные матрицы 4.2.9. Симметричный выбор ведущего элемента 4.3 Ленточные системы 141 Во многих приложениях, в которых присутствуют линейные системы, матрицы коэффициетов является ленточной. В этом случае уравнения всегда могут быть упорядочены так, что каждое неизвестное x(i) появляется только в нескольких уравнениях, соседних с i-м уравненинем. Формально мы говорим, что матрица A = (a(ij)) имеет верхнюю ширину ленты q, если a(ij) = 0 для всех j > i + q, и нижнюю ширину ленты p, если a(ij) = 0 для всех i > j + p. При решении ленточных систем можно добиться реальной экономии, так как треугольные множители LU, GG(T), LDM(T) и т.д., также являются ленточными. Прежде чем продолжить, мы советуем читателю просмотреть 1.2, где обсуждаются некоторые аспекты преобразований с ленточными матрицами. 4.3.1. Ленточное LU-разложение 4.3.2. Решение треугольных ленточных систем 4.3.3. Структура данных летночной матрицы 4.3.4. Ленточное исключение Гаусса с выбором ведущего элемента 4.3.5. LU-разложение матрицы Хессенберга 4.3.6. Ленточный метод Холецкого 4.3.7. Решение трехдиагональных систем 4.3.8. Вопросы векторизации 4.4 Симметричные неопределенные системы 150 Симметричная матрица, квадратичная форма которой x(T)Ax принимает и положительные, и отрицательные значения, называется неопределенной. Хотя неопределенная матрица может иметь LDL(T)-разложение, элементы множителей могут быть произвольной величины: [e 1] [1 0][e 0][1 0]T = [1 0] [1/e 1][0 -1/e][1/e 1] . Конечно, можно призвать на помощь некоторые стратегии с выбором из 3.4. Однако они теряют симметрию, и вместе с этим пропадает шанс для построения метода решения неопределенных систем "со скоростью Холецкого". Надо применить, согласно нашим рекомендациям из 4.2.9, симметричный выбор ведущего элемента, т.е. преобразование данных вида A <- PAP(T). К несчастью, симметричный выбор не всегда устойчив при вычислении LDL(T)-разложения. Если e(1) и e(2) малы, то независимо от матрицы P матрица [e(1) 1] A(~) = P [1 e(2)] имеет маленькие элементы на диагонали и в процессе разложения появляются большие числа. При симметричном выборе ведущие элементы всегда выбираются из диагонали и получаются ужасные результаты, если эти числа малы относительно тех внедиагональных элементов, которые обнуляются. Таким образом, LDL(T)-разложение не может быть рекомендовано как надежный подход к решению симметричной неопределенной системы. По-видимому, это действительно очень сложная задача - включить внедиагональные элементы в процессе выбора ведущего элемента и при этом сохранить симметрию. В данном разделе мы обсудим два способа, как это можно сделать. Первый метода предложен Aasen (1971), он вычисляет разложение PAP(T) = LTL(T), где L = ( l(ij) )-нижняя унитреугольная матрица, а T-трехдиагональная. Матрица P является перестановочной и выбираетс так, что |l(ij)| <= 1. Метода с диагональным выбором вычисляет матрицу P так, что PAP(T) = LDL(T), где матрица D состоит из блоков размеров 1 x 1 и 2 x 2. И снова матрица P выбирается так, что элементы нижней унитреугольной матрицы L удовлетворяют условию |l(ij)| <= 1. Оба разложения тебуют n**3/3 флопов и после вычисления могут использоваться для решения Ax = b с объемом работы O(n**2): PAP(T) = LTL(T), Lz = Pb, Tw = z, L(T)y = w, x = Py -> Ax = b, PAP(T) = LDL(T), Lz = Pb, Dw = z, L(T)y = w, x = Py -> Ax = b. В этих процедурах имеется только одна новая тема, которую имеет смысл обсудить - это решение системы Tw = z и Dw = z. В методе Аазена симметричная неопределенная трехдиагональнам система решается за O(n) операций при помощи исключения Гаусса с выбором. Отметим, что на этом уровне вычислений не будет серьезной платы за невнимание к симметрии, поскольку общий объем вычислений O(n**3). В подходе с диагональным выбором система Dw = z составлена из множества симметричных неопределенных систем порядка 1 x 1 и 2 x 2. Задача порядка 2 x 2 могут быть рашены посредством исключения Гаусса с выбором. Опять же не будет вреда от пренебрежения симметрией на этой фазе вычислений сложности O(n). Таким образом, центральный вопрос этого раздела - это эффективное вычисление разложений (4.4.1) и (4.4.2). 4.4.1. Алгоритм Парлетта-Рейда 4.4.2. Метод Аазена 4.4.3. Выбор ведущего элеметна в методе Аазена 4.4.4. Методы с диагональным выбором 4.4.5. Устойчивость и эффективность 4.4.6. Сравнение метода Аазена с диагональным выбором 159 (80) 4.5 Блочные трехдиагональные системы 159 Во многих областях приложения возникают матрицы, имеющие блочную структуру, которую можно эффективно использовать. Например, в задачах оптимизации с ограничениями часто требуется решать линейные системы вида [A B][y] [c] [ ][ ] = [ ], [B(T) 0][z] [d] где A - симметричная положительно определенная матрица, а B-матрица полного столбцового ранга. В данной ситуации имеет смысл использовать структуру системы, а не трактовать (4.5.1) как еще одну симметричную неопределенную систему. Для ознакомления с деталями см. Heath (1978). Чтобы изучить проблему, как именно можно использовать блочную структуру, выберем для рассмотрения решение блочных трехдиагональных систем вида [D1 F1 .. 0 ][x1] [b1] [E1 D2 .. .. . ][x2] [b2] [. . ][. ] = [ .] [. F(n-1)][. ] [ .] [0 ....... E(n-1) D(n)][xn] [bn] Здесь мы считаем, что все блоки имеют размер q x q и что x(i) и b(i) принадлежат R(q). В этом разделе мы рассмотрим два подхода к этой задаче - блочное LU-разложение и уступающую ему схему, известную как циклическая редукция. Другие аспекты, связанные с блочными трехдиагональными системами, обсуждаются в 10.3.3. 4.5.1. Блочное LU-разложение 4.5.2. Блочное диагональное доминирование 4.5.3. Сравнение блочного и ленточного решений 4.5.4. Блочная циклическая редукция 4.6 Системы Вандермонда 166 Пусть вектор x(0:n) -> R(n+1). Матрица V -> R(n+1)x(n+1) вида [1 1 ... 1] [x0 x1 ... xn] [. . ] V = V(x0,......,xn) = [. . ] [. . ] [xn(0) xn(1) .... xn(n)] называется матрицей Вандермонда. В этом разделе мы покажем, как можно решить системы V(T)a = f = f(0:n) и Vz = b = b(0:n) за O(n**2) флопов. Для удобства нижний индекс элементов векторов и матрицы начинается с 0. 4.6.1. Интерполяционный многочлен V(T)a = f 4.6.2. Система Vz = b 4.6.3. Устойчивость 4.7 Тёплицевы системы 171 (86) Матрицы, элементы которых равны вдоль каждой диагонали, возникают во многих приложениях и называются тёплицевыми матрицами. Формально матрица T -> R(n,n) является тёплицевой, если существуют скаляры r(n+1),....,r(),r(n-1), такие что a(ij) = r(j-i) для всех i и j. Таким образом, матрица [r0 r1 r2 r3] [r-1 r0 r1 r2] [r-2 r-1 r0 r1] [r-3 r-2 r-1 r0] является тёплицевой. Тёплицевы матрицы относятся к обширному классу персимметричных матриц. Говорят, что матрица B -> R(n,n) является персимметричной, если она симметрична относительно северо-юго-западной диагонали, т.е. b(ij) = b(n-j+1,n-i+1) для всех i и j. Это эквивалентно требованию B = E B(T) E, где E = [e(n),.....,e(1)] = I(n)(:,n:-1:1) является n x n-матрицей перестановок. Легко проверить, что (а) тёплицевы матрицы - персимметричные и (b) обратная к невырожденной тёплицевой матрице является персимметричной. В этом разделе мы покажем, как правильное использование свойства (b) поможет нам решать тёплицевы системы за время O(n**2). Обсуждение ограничивается важным случаем, когда матрица T симметричная и положительно определенная. 4.7.1. Три задачи Допустим, что у нас имеются скаляры r1,.....,rn, такие, что для k = 1:n матрицы [1 r1 ... r(k-2) r(k-1)] [r1 1 ... r(k-2)] [ . . . . ] [ . . . . ] T(k) = [ . . . . ] [ . . . . ] [r(k-2) . . r1] [r(k-1) r(k-2) ... r1 1] положительно определенные. (Мы не ограничиваем общности нормализацией диагонали). В этом разделе описаны три алгоритма: Алгоритм Дурбина для задачи Юла-Уолекра T(n)y = -(r1,....r(n))(T). Алгоритм Левинсона для задачи с произвольной правой частью T(n)x= b. Алгоритм Тренча для вычисления В = 1/T(n) При выводе этих методов обозначи через E(k) перестановочную k x k-матрицу, Е(k) = I(k)(:,k:-1:1). 4.7.2. Решение уравнений Юла-Уолкера Мы начинаем с представления алгритма Дурбина для уравнений Юла-Уолкера, которые возникают в связи с линейными задачами предсказания. Пусть для некоторое k, которое удовлетворяет условию 1 <= k <= n - 1, мы решили систему Юла-Уолкера T(k)y = -r = -(r1,...,r(k))(T) порядка k. Покажем, как теперь можно решить систему порядка k + 1 173 (87) [T(k) E(k)r][z] [r ] = [r(T)E(k) 1][a] [r(k + 1)] за O(k) флопов. Сначала заметим, что z = (-r - aE(k)r)/T(k) = y - aE(k)r/T(k) и a = -r(k + 1) - r(T)E(k)z Так как T(k)(-1) персимметричная, T(k)(-1)E(k) = E(k)T(k)(-1), и поэтому z = y - aE(k)T(k)(-1)r = y + aE(k)y. Подсталяя это в верхнее выражение для a, находим a = -r(k + 1) - r(T)E(k)(y + aE(k)y) = -(r(k + 1) + r(T)E(k)y)/(1 + r(T)y). Знаменатель является положительным, поскольку T(k + 1) полжительно определена и, кроме того, [I E(k)y] [T(k) E(k)r][I E(k)y] [T(k) 0 ] (T) = [0 1 ] [r(T)E(k) 1][0 1 ] [0 1 + r(T)y]. Мы продемонстрировали k-й шаг алгоритма, предложенного Дурбиным (1960). Он сотоит в решении Юла-Уолкера T(k)y(**k) = -r(**k) = -(r1,.....,r(k))**T для k = 1:n и может быть записан следующим образом: y**1 = -r(1) for k = 1:n - 1 b(k) = 1 + [r**k]**T y**k a(k) = - (r(k+1) + r(**k)**T E(k)y**k/b(k) z(k) = y**k + a(k)E(k) y**k [z**k] y**(k + 1) = [a(k)] end Как теперь нами установлено, этот алгоритм требует 3n**2 флопов для произвольного вектора y = y**(n). Возможно, однако, уменьшить общий объем работы, используя некоторые приведенные выше выражения: b(k) = 1 + [r**(k)]**T y**(k) = = 1 + [r**(k-1)**T r(k)][y**(k-1) + a(k-1)E(k-1)y**(k-1)] = = (1 + [r**(k-1)]**T y(k-1)) + a(k-1)( [r**(k-1)]**T E(k-1)y**(k-1) + r(k) ) = = b(k-1) + a(k-1)(-b(k-1)a(k-1)) = = (1 - a(k-1)**2)b(k-1). Применяя эту рекурсию, мы получаем следующий алгоритм: Алгоритм 4.7.1 (Дурбин). Даны вещественные числа 1 = r0, r1, ....., rn, такие, что матрица T = ( r(|i-j|) ) -> R**(n x n) является положительно определенной, следующий алгоритм вычисляет вектор y -> R**n, такой, что Ty = -(r1,...,rn)**T. y(1) = -r(1); b = 1; a = -r(1) for k = 1:n - 1 b = (1 -a**2)b a = -(r(k+1) + r(k:-1:1))**T y(1:k))/b for i = 1:k z(i) = y(i) + ay(k + 1 - i) end y(1:k) = z(a:k); y(k+1) = a end Алгоритм требует 2n**2 флопов. Мы включили вспомогательный вектор z для прозрачности алгоритма, но он может быть убран. Пример 4.7.1. Пусть мы хотим решать систему Юла-Уолекра [1 0.5 0.2][y1] [0.5] [0.5 1 0.5][y2]=-[0.2], [0.2 0.5 1][y3] [0.1] используя алгоритм 4.7.1. После первого происхождения через цикл мы получим [-8/15] a = 1/15, b = 3/4, y = [ 1/15] Потом вычислим величины b = (1 -a**2)b = 56/75, a = -(r3 + r2y1 + r1y2)/b = -1/28, z1 = y1 + ay2 = -225/420, z2 = y2 + ay1 = -36/420 дающие окончательное решение y = (-75,12,-5)**T/140. 174 88 4.7.3. Задача с произвольной правой частью Если еще немножко потрудиться, можно решить симметричную положительно определенную тёплцеву систему, которая имеет произвольную правую часть. Предположим, что ма решили систему T(k)x = b = (b1,....,bk)**T (4.7.2) для некоторого k, удовлетворяющего условию 1 <= k < n, и что сейчас мы хотим решить систему [T(k) E(k)z][ v ] [b ] = (4.7.3) [r**T*E(k) 1][ m ] [b(k+1)] Здесь, как и выше, r = (r1,....,n)**T. Допустим также, что решение систему Юла-Уолкера k-го порядка T(k)y = -r уже получено. Так как v = T(k)**(-1)(b - mE(k)r) = x + mE(k)y, то отсюда следует, что m = b(k+1) - r**T*E(k)v = b(k+1) - r**T*E(k)x - mr**Ty =(b(k+1) - r**T*E*(k)*x)/(1 + r**T*y). Следовательно, мы можем сделать переход от (4.7.2) к (4.7.3) за O(k) флопов. В целом мы можем эффективно решить систему T(n)x = b, решая системы T(k)x**k = b**k = (b1,....,bk)**T и T(k)y**k = -r**k = (r1,....,rk)**T "параллельно" для k = 1:n. Это составляет содержание следующего алгоритма. Алгоритм 4.7.2.(Левинсон). Даны вектор b -> R(n) и вещественные числа 1 = r0, r1,....., rn, такие, что матрица T( r(|i-j|) ) -> R(n,n) является положительно определенной; следующий алгоритм вычисляет вектор x -> R(n), такой, что Tx = b. y(1) = -r(1); x(1) = b(1); b = 1; a = -r(1) for k = 1:n - 1 b = (1-a)**2*b; m = (b(k + 1) - r(1:k)**T*x(k:-1:1))/b v(1:k) = x(1:k) + my(k:-1:1) x(1:k) = v(1:k); x(k:-1:1) = m if k < n - 1 a = -(r(k+1) + r(1:k)**T*y(k:-1:1))/b z(1:k) = y(1:k) + ay(k:-1:1) y(1:k) = z(1:k); y(k+1) = a end end Алгоритм требует 4n**2 флопов. Векторы z и v введены для ясности, и без них можно обойтись. Пример 4.7.2. Пусть мы хотим решить симметричную положительно определенную тёплицеву систему [ 1 0.5 0.2][x1] [ 4] [0.5 1 0.5][x2]= [-1], [0.2 0.5 1 ][x3] [ 3] используя приведенный выше алгритм. После первого прохождения через цикл мы получим [-8/15] [ 6] a = 1/15, b = 3/4, y =[ ], x = [ ]. [ 1/15] [-4] Потом вычислим величины b = (1 - a**2)b = 45/75, m = (b3 - r1x2 - r2x1)/b = 285/56, v1 = x1 + my2 = 355/56, v2 = x2 + my1 = -376/56, дающие окончательное решение x = (355, -376, 285 )**T/56. 4.7.4. Вычисление обратной матрицы Одно из наиболее удивительных свойств симметричной положительно определенной тёплицевай матрицы T(n) состоит в том, что обратная к ней может быть вычислена полностью за O(n**2) флопов. Чтобы получить алгоритм, реализующий эту процедуру, разобьем T(n)**-1 следующим образом: [A Er ] [B v] T(n)**-1 = [ ] = [ ] [r**T*E 1 ] [v**T gamma] где A = T(n-1), E = E(n-1) и r = (r1,...,r(n-1))**T. Из уравнения [A Er][ v] [0] [ ][ ]= [ ] [r**T*E 1][gamma] [1] следует, что Av = -gammaEr = -gammaE(r1,....,r(n-1))**T и gamma = 1 - r**T*E*v. Если вектор y разрешает систему Юла-Уолкура (n - 1) порядка, Ay = -r, тогда из этих выражений вытекает, что gamma = 1/(1 + r**T*y), v = gamma*E*y. Таким образом, последние строку и столбец матрицы T(n)**-1 получить легко. Нам осталось развить формулы для получения элементов подматрицы B из (4.7.4). Так как AB + E*r*v**T = I(n - 1) , отсюда следует, что B = A**-1 - (A**-1*E*r)v**T = A**-1 + v*v**T/gamma. И так как A = T(n - 1) является невырожденной и тёплицевой, обратная к ней персимметрична. Таким образом, b(ij) = (A**-1)(ij) + v(i)v(j)/gamma = = (A**-1)(n-j,n-i) - v(n-j)v(n-i)/gamma + v(i)v(j)/gamma = = b(n-j,n-i) + (v(i)v(j) - v(n-j)v(n-i))/gamma Это показывает, что, хотя B не является персимметричной, мы можем легко вычислить элемент b(ij), используя его отражение относительно северо-восточной-юго-западной диагонали. Это в сочетании с фактом, что A**-1 персимметричная, приводит нас к определению B от "краев" к "середине". Поскольку порядок операций правильно описать затруднительно, мы сначала покажем наглядно формальные детали алгоритма. С этой целью допустим, что известны последние строка и столбец матрицы A**-1: [u u u u u k] [u u u u u k] [u u u u u k] A**-1 = [u u u u u k] [u u u u u k] [k k k k k k] Здесь u и k обозначают неизвестные и известные элементы соотвественно и n = 6. Поочередно используя персимметрию матрицы A**-1 и рекурсию (4.7.5.), мы можем вычислить B следующим образом: [k k k k k k] [k k k k k k] [k u u u u k] [k u u u k k] [k u u u u k] [k u u u k k] персимм. [k u u u u k] (4.7.5) [k u u u k k] ---> [k u u u u k] ----> [k k k k k k] [k k k k k k] [k k k k k k] [k k k k k k] [k k k k k k] [k k k k k k] [k k k k k k] [k k u u k k] [k k u k k k] персимм. [k k u u k k] (4.7.5) [k k k k k k] ---> [k k k k k k] ----> [k k k k k k] [k k k k k k] [k k k k k k] [k k k k k k] [k k k k k k] [k k k k k k] персимм. [k k k k k k] ---> [k k k k k k] [k k k k k k] Конечно, когда вычисляемая матрица одновременно и симметричная и персимметричная, такая, как A**-1, то необходимо вычислять только "верхний клин" матрицы, например x x x x x x x x x x x x (n = 6). С учетом последнего замечания, мы готовы представить алгоримть полностью. Алгоритм 4.7.3. (Тренч). Даны вещественные числа 1 = r0,r1,....,rn , такие, что матрица T = (r(|i-j|)) -> R(n,n) является положительно определенной, следующий алгоритм вычисляет матрицу B = T(n)**-1. Вычисляются только те b(ij), для которых i <= j и i + j <= n + 1. Используя алгоритм 4.7.1, решить систему T(n-1)y = -(r1,....,r(n-1))**T . gamma = 1/(r(1:n-1)**Ty(n:-1:1)) v(1:n-1) = gamma y(n-1:-1:1)**T B(1,1) = gamma B(1,2:n) = v(n-1:-1:1)**T for i = 2:floor((n-1)/2) + 1 for j = i:n - i + 1 B(i,j) = B(i-1,j+1) + (v(n+1-j)v(n+1-i) - v(i-1)v(j - 1))/gamma end end Алгоритм требует 13n**2/4 флопов. Пример . Если приведенный выше алгоритм применяется для вычисления обратной матрицы B полжительно определенной тёплицевой матрицы [ 1 0.5 0.2] [0.5 1 0.5] [0.2 0.5 1 ] то мы получим gamma = 75/56, b11 = 75/56, b12 = -5/7, b13 = 5/56 и b22 = 12/7. 178 90 4.7.5. Вопросы устойчивости Анализ ошибок упомянутых выше алгоритмов был провден Cybenko (1978), и мы вкратце опишем некоторые из его результатов. Оказывается, ключевыми величинами являются а(k) из (4.7.1). В точной арифметике эти скаляры удовлетвояют соотношению |a(k)| < 1 и могут быть использованы для оценки ||T**-1||(1): max(1/П(j=1,n-1)(1 - a(j)**2), 1/П(j=1,n-1)(1 - a(j))) <= ||T(n)**-1|| <= П(j=1,n-1)((1 + |a(j)| )/(1 - |a(j)| )) Более того, решение системы Юла-Уолкера T(n)y = -r(1:n) удовлетворяет соотношению ||y||(1) = ( П(k=1,n-1)(1 + a(k))) - 1 в том случае, если все a(k) неотрицательные. Теперь если x(^) - это вычисленное алгритмом Дурбина решение уравнений Юла-Уолкера, то r(p) = T(n)x(^) + r может быть оценено следующим образом: ||r(p)|| = u П(k=1,n)(1 + |a(k)(^)|), где a(k)(^) - вычисленный аналог a(k). С помощью сравнения, так как каждый |r(i)| ограничен единицией, получаем, что ||r(c)|| = u ||y||(1), где r(c) -невязка, соответсвующая вычисленному решению, полученному алогритмом Холецкого. Заметим, что невязки могут быть сравнены по величине в том случае, если справедлива формула (4.7.7). Эксперимента дает основания предположить, что это имеет место, даже если некоторые из a(k) отрицательные. Аналогичныц объяснения применимы к численному поведению алгоритма Левинсона. Для метода Тренча можно показать, что матрица B - вычисленный аналог обратной к T(n)**-1, удовлетворяет условию ||T(n)**-1 - B(^)||/||T(n)(-1)|| = u П()((1 + |a(k)(^)|)/(1 - |a(k)(^)|)). В свете формулы (4.7.7) видно, что правая часть является приближенной верхней оценкой для u||T(n)**-1|| , которая аппроксимирует величину относительной ошибки, когда T(n)**-1 вычислена при помощи разложения Холецкого. Замечания и литература к 4.7. Всякий кто отважится проникнуть в литературу по быстрым тёплицевым методам, должен прежде всего прочитать Bunch J.R.(1985). "Stability of Methods for Solvinf Toeplitz systems of Equations", SIAM J.Sci.Stat.Comp. 6, 349-364. для выяснения вопросов устойчивости. Как и вообще, среди быстрых алгоритмов есть много неустойчивых тёплицевых методов и необходимо проявлять осторожность. Первоначальными работами по трем методам, описанным в этом разделе, являются Durbin J.(1960)."The Fitting of Time Series Models", Rev.Inst.Int.Stat.28,233-243. Levinson N.(1947)."The Weiner RMS Error Criterion in Filter Design and Prediction", J.Math.Phys. 25,261-278. Trench W.F.(1964)."An Algorithm for the Inversion of Finite Toeplitz Matrices", J.SIAM 12, 515-522. Более детальное описание несимметричного алгоритма Тренча даетс в работе Zohar S.(1969)."Toeplitz Matrix Invertion: The Algorithm of W.F. Trench",J.ACM 16, 592-601. Другие ссылки относительно обращений тёплицевой матрицы содержит работы Trench W.F.(1974)."Invertion of Toeplitz Band Matrices",Math.Comp.28, 1089-1095. Watson G.A.(1973)."An Algorithm for the Invertion of Block Matrices of Toeplitz Form", J.ACM 20,409-415. Оценки ошибок, которые мы приводили, взятыRRизализа ошибок округления, имеюдегося в работах Cybenko G.(1978)."Error Analysis of Some Signal Processing Algorithms", Ph.D.thesis, Princeton Untversity. Cybenko G.(1980)."The Numerical Stabilbity of the Levinson-Durbin Algorithm fo Toeplitz Systems of Equations", SIAM J.Sci. and Stat.Comp. 1, 303-319. Также известны методы треугольного разложения тёплицевых систем сложности O(n**3). Phillips J.L.(1971)."The Trinangulat Decomposition of Hankel Matrices", Math.Comp. 25, 599-602. Rissansen J.(1973)."Algorithms for Triangular Decomposition of Block Hankel and and Toeplitz Matrices with Application for Factoring Positive Matrix Polynomials", Math.Comp. 27, 147-154. Упомянем следующие работы, описывающие некоторые важные приложения тёплицевых матриц: Makhoul J.(1975)."Linear Prediction: A Tutoial Review", Proc. IEEE 63(4), 561-580. Markel J. and Gray A.(1976).Linear Prediction of Speech, Springer-Vetag, Berlin and New York. Oppenheim A.V.(1978). Application of Digital Signal Processing, Prentice-Hall, Englwood Cliffs. Глава 5. Ортогонализация и метод наименьших квадратов 181 В этой главе речь пойдет главным образом о решении переопределенных систем уравнений методом наименьших квадратов, т.е. о минимизации функционала ||Ax - b||2, где A -> R(m,n), m >= n, b -> R(m). Наиболее надежные методы решения этой задачи связаны с приведением матрицы A к разнообразным каноническим формам путем ортоганальных преобразований. Центральную часть при этом играют отражения Хаусхолдера и вращения Гивенса, с рассмотрения которых мы и начинаем эту главу. В 5.2. обсуждаются вычисление разложения A = QR, где Q - ортогональная, R - верхняя треугольная матрица. Оно сводится к нахождению отртонормального базиса в range(A). В 5.3 показано, что QR-разложение может быть применено для решения задачи мнимальных квадратов с матрицей полного ранга. Этот метод мы затем сравниваем с методом нормальных уравнений, предварительно развив соотвествующую теорию возмущений. В 5.4. и 5.5 мы рассматриваем методы, пригодные в сложной ситуации, когда матрица A неполноранговая (или почти неполноранговая). Здесь представлены QR-разложение с выбором ведущего элемента по столбцу и SVD. В 5.6 мы обсуждаем некоторые шаги, которые могут быть предприняты для улучшения качества вычисленного решения задачи минимальных квадратов. 5.7 содержит несколько замечаний о недоопределенных системах. 5.1 Матрицы Хаусхолдера и Гивенса 181 Напомним, что -матрица называется ортогональныой, если . Ортогональные матрицы играют важную роль при решении задач на собственные значения и задач минимальных квадратов. В этом разделе мы вводим ключевые для эффективных вычислений с ортогональными матрицам преобразования: отражения Хаусхолдера и вращения Гивенса. 5.1.1. Двумерный случай Полезно посмотреть, какой геометрический смысл имеют вращения и отражения в случае n = 2. Ортогональная 2 x 2-матрица Q называется матрицей вращения, если [cos(fi) sin(fi)] Q = [ ], [-sin(fi) cos(fi)] Если y = Q**T*x , то y получается поворотом x на угол fi против часовой стрелки. Ортогональная 2 x 2-матрица называется матрицей отражения, если [cos(fi) sin(fi)] Q = [ ]. [sin(fi) -cos(fi)] Если y = Q**T, x = Qx, то y получается отражением x относительно оси, определяемой как [cos(fi)/2] S = span [ ] [sin(fi)/2] Вращения и отражения привлекательны с вычислительной точки зрения, так как они легко строятся и могут быть использованы для обнуления отдельных элементов векторов при надлежащем выборе угла поворота или плоскости отражения. 182 (92) Пример 5.1.1. Пусть x = (1, root(3))**T. Положим [cos(-60 гр) sin(-60 гр)] [ 1/2 - root(3)/2] Q = [ ] = [ ] [-sin(-60 гр) cos(-60 гр)] [root(3)/2 1/2 ] Тогда Q**T*x = (2,0)**T. Таким образом, поворот на -60 градусов обнуляет вторую компоненту X. Если [cos(30 гр) sin(30 гр)] [root(3)/2 1/2 ] Q = [ ] = [ ] [sin(30 гр) -cos(30 гр)] [ 1/2 -root(3)/2] то Q**T*x = (2,0)**T. Таким образом, вторую компоненту x можно обнулить, отражая x относительно прямой с наклоном 30 градусов. 5.1.2. Отражение Хаусхолдера 5.1.3. Вычисление вектора Хаусхолдера 5.1.4. Умножение на матрицы Хаусхолдера 5.1.5. Ошибки округления 5.1.6. Факторизованное представление 5.1.7. Блочное представление 5.1.8. Вращение Гивенса 5.1.9. Умножения на матрицы Гивенса 5.1.10. Ошибки округления 5.1.11. Представление произведений матриц Гивенса 5.1.12. Распространение ошибок 5.1.13. Быстрые вращения 5.2 QR-разложение 195 В этом разделе мы покажем, как применить преобразования Гивенса и Хаусхолдера для получения различных разложений матрицы. Мы начинаем с QR-разложения, которое для m x n-матрицы A выглядит так: A = QR, где Q -> R(m,m) - ортогональная, R -> R(m x n) - верхняя треугольная матрица. В этом разделе мы предполагаем m >= n. Если A имеет полный столбцовый ранг, то первые n столбцоа матрицы Q образуют ортонормированный базиc подпространства range(A). Тем самы QR-разложение дает один из способов получения ортонормированного базиса для набора векторов. Существует несколько путей вычисления QR-разложения. Мы изложим методы, основанные на преобразованиях Хаусхолдера и Гивенса, а также быстрых вращениях и блочных отражениях. Кроме того, мы обсуждаем процесс ортогонализации Грамма-Шмидта и его модификацию, отличающуюся большей численной устойчивостью. 5.2.1. QR-разложение: преобразования Хаусхолдера 5.2.2. QR-разложение: метод блочных отражений 5.2.3. QR-разложение: преобразования Гивенса 5.2.4. QR-разложение хессенберговской матрицы с использоанием вращений 5.2.5. QR-разложение: быстрые вращения 5.2.6. Свойства QR-разложения 5.2.7. Классический метод Грама-Шмидта 5.2.8. Модифицированный метод Грама-Шмидта 5.2.9. Вычислительные затраты и точность 5.3 Задача наименьших квадратов: случай полного ранга 205 Рассмотрим задачу нахождения вектора x, удовлетворяющего уравнению Ax = b, где матрица данных A -> R(m,n) и вектор наблюдений b -> R(m) заданы, причем m >= n. Когда уравнений больше, чем неизвестных, мы говорим, что система Ax = b переопределена. Обычно переопределенная система не имеет точного рашения, поскольку b должен принадлежать range(A), являющемуся собственным подпространством R(m). Чтобы все же как-то решить задачу, попытаемся минимизировать ||Ax - b||(p) для некоторого подходящего p. Оптимальные решения будут различным в разных нормах. Например, если A = [1 1 1]**T, b = (b1 b2 b3)**T, где b1 >= b2 >= b3 >= 0, то можно проверить, что p = 1 -> x(opt) = b(2), p = 2 -> x(opt) = (b1 + b2 + b3)/3, p = 00 -> x(opt) = (b1 + b3)/2. Минимизация в 1-норме и в 00-норме затруднена тем, что функция f(x) = ||Ax -b||(p) не дифференцируема для этих значений p. Несмотря на это, в этой области был достигнут существенный прогресс, и для миниминзации в 1-норме и в 00-норме имеется несколько хороших методик. См. Bartels, Conn, Sinclair (1978), Bartels, Conn, Charalambous (1978). В отличие от минимизации в произвольной p-норме, задача наименьших квадратов (задача LS) min||Ax - b||(2) (5.3.1) лучше поддается рашению по двум причинам: w(x) = 1/2||Ax - b||(2)**2 - дифференцируемая функция x, поэтому значение x в точке минимума удовлетворяет градиетному уравнению grad(w(x)) = 0. Это уравнение оказывается легко строящейся симметричной линейной системой, положительно определенной, если A имеет полный столбцовый ранг. 2-норма сохраняется ортогональными преобразованиями. Это означает, что мы имеем право искать такую ортогональную матрицу Q, что эквивалентная задача минимизации ||(Q**T*A)x - (Q**T*b)||(2) решается "легко". В этом разделе мы проследим эти два подхода к решению для случая, когда A имеет полный столбцовый ранг. Будут подробно рассмотрены и сопоставлены методы, основанные на нормальных уравнениях и QR-разложении. 5.3.1. Следствие полноты ранга 5.3.2. Обусловленность прямоугольных матриц 5.3.3. Метод нормальных уравнений 5.3.4. Решение задач LS через QR-разложение 5.3.5. Срыв в случае "почти неполноранговости" 5.3.6. О методе MGS 5.3.7. Решение задачи LS методом быстрых вращений 5.3.8. Чувствиетльность задачи LS к возмущениям 5.3.9. Сравнение метода нормальных уравнений и QR-разложения 5.4 Другие ортогональные разложения 215 Если матрица A неполного ранга, то QR-разложение необязательно дает базис подпространства range(A). Эта проблема может быть решена, если перед QR-разложением переставить столбцы матрицы A, т.е. AП = QR, где П - матрица перестановок. Данные в матрице A могут быть сконцентрированы еще больше, если разрешить умножение справа на подходящую ортогональную матрицу Z:Q**T*AZ = T. Для выбора Q и Z имеются интересные возможности, которые, наряду с QR-разложением с выбором ведущего столбца, будут рассмотрены в этом разделе. 5.4.1. Случай неполного ранга: QR с выбором ведущего столбца 5.4.2. Полные ортогональные разложения 5.4.3. Двухдиагонализация 5.4.4. R-двухдиагонализация 5.4.5. Сингулярное разложение и его вычисление 5.5 Задача LS неполного ранга 221 Если матрица A неполного ранга, то задача LS имеет бесконечное множество решений, и мы вынуждены прибегать к специальным приемам, направленным на решение трудной задачи численного определения определения ранга. После некоторых предварительных замечаний, связанных с SVD, мы показываем, как можно использовать QR-разложение с выбором ведущего столбца для нахождения решения x(b) с тем свойством, что A*x(b) есть линейная комбинация r = rank(A) столбцов. Затем мы обсуждаем решение с минимальной 2-нормой, которое может быть получено через SDV. 5.5.1. Решение с минимальной нормой 5.5.2. Полное ортогональное разложение и x(ls) 5.5.3. Сингулярное разложение и задача LS 5.5.4. Псевдообратная матрица 5.5.5. Некоторые вопросы чувствительности к возмущениям 5.5.6. QR с выбором ведущего столбца и основные решения 5.5.7. Численное нахождение ранга с помощью aП = QR 5.5.8. Численный ранг и SVD 5.5.9. Некоторые сравнения 5.6 Взвешивание и итерационное уточнение 230 В 3.5. в контексте квадратных линейных систем были введены понятия масштабирования и итерационного уточнения. Здесь мы обобщаем эти идеи на случай задачи наименьших квадратов. 5.6.1. Взвешивание по столбцам 5.6.2. Взвешивание по строкам 5.6.3. Обобщенные наименьшие квадраты 5.6.4. Итерационное уточнение 5.7 Квадратные и недоопределенные системы 234 Развитые в этой главе метода ортогонализации могут быть применены к решешию квадратных систем, а также систем, в которых уравнений меньше, чем неизвестынх. В этом коротком разделе мы обсудим некоторые из возможнотей. 5.7.1. Использование QR и SVD систем для решения квадратных систем 5.7.2. Недоопределенные сисмемы 5.7.3. Возмущенные недоопределенные системы Глава 6. Параллельные матричные вычисления 238 Параллельные вычисления с матрицами - это область интенсивных исследований. Поскольку данное научное направление зародилось совсем недавно, литературы здесь предоставлено и в ней преобладает "изучение отдельных случаев". Трудности машинной зависимости мы обходим, занимаясь проработкой алгоритмов на достаточно высоком уровне. Рассматриваются прадигмы распределенной и общей памяти. Технические детали сведены к минимуму, и мы лишь слегка касаемся нескольких важных языковых и системных проблем, чтобы сосредоточиться на главных вопросах. Как следствие, в этой главе мало конкретных рецептов - здесь вы не найдете четких рекомендаций о том, как решать задачу X на параллельной машине Y. В 6.1 и 6.2 ключевые идеи вводятся на примере операции gaxpy. В 6.3 обсуждается параллельное умножение матриц для систем как с распределенной, так и с общей памятью. Затем в 6.4, 6.5 и 6.6 рассматривается паралльельное вычисление различных матричных разложений. Мы делаем акцент на разложении Холецкого, но всюду даются указания насчет LU- и QR-разложений. 6.1 Операции на распределенной памяти 238 В этом разделе разрабатываются средства записи алгоритмов для распределенной памяти. Чтобы проиллюстрировать их и выявить ключевые положения, рассматриыается операция gaxpy: z = y + Ax, A -> R(n,n), x,y,z -> R(n) Более сложные вычислительные процессы будут рассматрены в 6.3 - 6.5. 6.1.1. Системы с распределенной памятью 6.1.2. Сети процессоров 6.1.3. Имена для соседей 6.1.4. Инициализация и окончание 6.1.5. Связь между процессорами 6.1.6. Некоторые распределенные структуры данных 6.1.7. Систолическая модель 6.1.8. Атрибуты параллельного алгоритма 6.1.9. Равномерная загруженность 6.1.10. Стоимость обмена информацией 6.1.11. Эффективность и ускорение 6.1.12. Дробление вычислений 6.1.13. Модель передачи сообщений 6.1.14. Дальнейшие примеры передачи сообщений 6.2 Операции на общей памяти 252 Теперь мы обсудим операцию gaxpy на многопроцессорных сисмтемах с общей памятью. В такой системе каждый процессор имеет доступ к общей, разделенной с другими памяти, как показано на следующем рисунке: Рис. 6.2.1. 4-проыессорная система с общей памятью Связь между процессорами достигается с помощью считывания и записи общих переменных, размещаемых в общей памяти. Каждый процессор выполняет свою локальную программу и имеет свою локальную память. В ходе вычислений данные текут туда и обратно из общей памяти. Например, для операции gaxpy x, y и A могут располагаться в общей памяти. Процессор, участвующий в обработке, может: a) скопировать x, i-ю блочную строку A(i) и i-й подвектор y(i) из общей памяти в свою память; b) вычислить z(i) = y(i) + A(i)x; с) скопировать z(i) из локальной памяти в общую память. В модифицированной форме здесь присутствует все то, чем мы интересовались в случае вычислений на распределенной памяти. Как обсуждалось в предыдущем праграфе, вся процедура должна обеспечивать равномерную загруженность. Вычисления следует организовывать таким образом, чтобы отдельным прдцесссорам приходилось как можно меньше простаивать в ожидании какой-нибудь полезной рабоыт. Потоки данных между общей памятью и локальными должны управляться с особой аккуратностью, так как они, по нашему мнению, являются значительным накладным расходом. (Это соотвествует межпроцессорным коммуникациям в случае распределенной памяти.) Природа физического соединения между процессорами и общей памятью очень важна и может влиять на разработку алгоритмов. Например, если в любой заданный момент времени доступ к общей памяти может получить только один процессор, то следует позаботиться о том, чтобы свести к минимуму возможность конфликтов, связанных с одновременным обращением к общей памяти. Большей частью мы стараемся не опускаться до столь подробной детализации, предпочитая рассматривать этот аспект системы как черный ящик, изображенный на рис. Все локальные вычисления подчинены традиционным требованиям эффективности. Если отдельные процессоры являются векторными/конвейерными процессорами, то, как обсуждалось в 1.4, могут быть важными таки вещи, как шаг и длина вектора, число обращений к локальной памяти повторного использования данных. Чтобы познакомится с матричными вычислениями на общей памяти, полезно еще раз рассмотреть параллельную реализацию для gaxpy. Сначала мы изложим такой алгоритм для gaxpy, в котором задача каждому процессору ставится заранее. Для менее регулярных вычислений необходимо назначить задачи процессорам по ходу вычислений. Динамическое назначение задач требует специальной синхронизации, и для этого мы вводим понятие монитора. Используя мониторы, мы может найти элегантные решения для широкого класса проблем планирования вычислений. 6.2.1. Статическое планирование операции gaxpy 6.2.2. Обмен данными с общей памяти 6.2.3. Равномерная загруженность 6.2.4. Синхронизация с помощью барьеров 6.2.5. Парадигма динамического резерва задач 6.2.6. Мониторы 6.2.7. Столбцовая реализация операции gaxpy 6.3 Параллельное умножение матриц 262 В этом параграфе рассматривается задача умножения матрицы на матрицу D = C + AB для многопроцессорных систем с распределенной и общей памятью. Рассматриваются процедуры для блочной операции gaxpy и для блочного скалярного произведения. 6.3.1. Процедуры для блочной операции 6.3.2. Проблема равномерной загруженности 6.3.3. Некоторые вопросы, связанные с дроблением вычислений 6.3.4. Систолическое умножение 3 x 3-матриц 6.3.5. Общий случай 6.3.6. Блочный аналог 6.3.7. Асинхронная тороидальная процедура 6.3.8. Использование блочного скалярного произведения в режиме резерва задач 6.4 Кольцевые процедуры разложения 273 В этом параграфе мы обсудим, как можно организовать различные матричные разложения на кольце. Представив параллельную реализацию основанного на gaxpy Холецкого, мы покажем затем, как работает та же самая логика планирования для разложений PA = LU и A = QR. В конце обсуждается параллельное решение треугольных систем. 6.4.1. Кольцевой алгоритм Холецкого: n = p 6.4.2. Кольцевой метод Холецкого: общий случай 6.4.3. Кольцевые процедуры для других разложений 6.4.4. Параллельное решение треугольных систем 6.5 Сеточные процедуры разложения 281 В этом параграфе обсуждаются различные алгоритмы для разложений Холецкого и QR. Для разложения Холецкого приводятся синхронные и асинхронные процедуры, позволяющие нам сравнить особенности рвссуждений, связанных с этим двумя способами вычислений на распределенных системах. Параграф заканчивается кратким обсуждением систолической сеточной процедуры для QR-разложения. 6.5.1. Сеточная систолическая процедура Холецкого 6.5.2. Асинхронная процедура Холецкого 6.5.3. Сеточная систолическая процедура QR-разложения 6.6 Методы разложения на общей памяти 290 В этом параграфе мы посмотрим на различные реализации разложения Холецкого на общей памяти. Каждая реализация обладает своим собственным обрзцом синхронизации, а вместе оин дают представление о том, как разложение типа Холецкого может организовываться на системе с общей памятью. Мы выбрали разложение Холецкого ради удобства изложения, но большей частью то, о чем мы говорим для него, переносится и на другие разложения из нашего репертуара. 6.6.1. Статическое планирование для Холецского с внешними произведениями 6.6.2. Статическое планирование для других разложений 6.6.3. Две параллельные реализации для Холецкого с gaxpy 6.6.4. Замечания об обменах с общей памятью 6.6.5. Блочный Холецкий на базе резерва задач Глава 7. Несимметричная проблема собственных значений 299 Рассмотрев линейные уравнения и метод наименьших квадратов, обратим внимание на третью важную область матричных вычислений - алгебраическую проблему собственных значений. В этой главе рассматривается несимметричная проблема, а более приятный симметричный случай - с следующей. Наша первоочередная задача - познакомится с разложениями Шура и Жордана и основными свойствами собственных значений и инвариантных подпространств. Противоположное поведение этих двух разложений положено в основу той части 7.2, в которой мы исследуем, как на собственные значения и инвариантные подпростраства влияет возмущение. Введены числа обусловленности, что позволяет получать оценку погрешности, которую можно ожидать из-за округления. Основной алгоритм данной главы - заслуженно известный QR-алгоритм. Эта процедура является наиболее сложным алгоритмом, рассмотренным в этой книге, и ее изложение занимает три раздела. Мы получаем базовые QR-итерации в 7.3 как естественное обобщение простейшего степенного метода. Последующие два раздела посвящены тому, чтобы сделать базовые итерации вычислительно осуществимыми. Это включает введение хессенбергова разложения в 7.4 и идею начальных сдвигов в 7.5. QR-алгоритм вычисляет вещественную форму Шура матрицы - каноническую форму, которая выявляет собственные значения, не не собственые векторы. Поэтому обычно необходимо выполнить дополнительные вычисления, если требуется информация, касающаяся инвариантных подпространств. В 7.6, который мог бы иметь подзаголово: "Что делать после того, как найдена вещественная форма Шура", мы обсуждаем различные методы вычисления инвариантных подпространств, которые позволяют довести QR-алгоритм до конца. Наконец, в последнем разделе мы рассматриваем обобщенную проблему собственных значений Ax = laymbdaBx и вариант QR-алгоритма, придуманный для ее решения. Этот алгоритм, называемый QZ-алгоритмом, подчеркивает важность ортогональных матриц в проблеме собственных значение - главную тему этой главы. Сейчас самое время сделать замечание о комплексной арифметике. В этой книге мы сосредоточиваем внимание на развитии вещественных алгоритмов для вещественных задач. Данная глава не представляет исключения, хотя даже вещественная несимметричная матрица может иметь комплексные собственные значения. Однако при получении практического вещественного QR-алгоритма и при математическом анализу самой проблемы собственных значений удобнее работать в поле комплексных чисел. Поэтому читатель обнаружит, что мы перешли к комплексным обозначениям в 7.1, 7.2 и 7.3. В этих разделах мы используем комплексное QR-разложение (A = QR, Q-унитарная, R-комплексная верхняя треугольные матрицы) и комплексные SVD (A = U Sigma V(H), U и V-унитарные, Sigma - вещественная диагонильная матрица). Получение этих комплексных разложений отличается от вещественного случая не так сильно, чтобы оправдать отдельное описание. 7.1 Свойства и разложения 300 В этом разделе мы даем обзор математических основ, неоходимых для построения и анализа обсуждаемых далее алгоритмов вычисления собственных значений. 7.1.1. Собственные значения и инвариантные подпространства 7.1.2. Основные унитарные разложения 7.1.3. Неунитарные преобразования 7.1.4. Некоторые замечания о неунитарном подобии 7.2 Теория возмущения 307 Никакие из разложений, описанных в предыдущем разделе, не могут быть вычислены точно из-за ошибок округления и из-за того, что алгоритмы вычисления собственных значений являются итерационными и должны быть остановлены после некоторого конечного числа шагов. Поэтому важно разработать удобную теорию возмущений, чтобы руководствоваться ею в последующих разделах. 7.2.1. Чувствительность собственного значения 7.2.2. Обусловленность простого собственного значения 7.2.3. Чувствительность кратных собственных значений 7.2.4. Чувствительность собственного вектора 7.2.5. Чувствительность инвариантного подпространства 7.3 Степенные итерации 316 Пусть заданы матрица A -> C(n,n) и унитарная матрица U(0) -> C(n,n). Предположим, что ортогонализация Хаусхолдера (алгоритм 5.2.1) может быть распространена на комплексные матрицы (это возможно), и рассмотрим следующий итерационный процесс T(0) = U(0)**H*A*U(0) for k = 1,2,.... T(k-1) = U(k)R(k) (QR-разложение) T(k) = R(k)U(k) end Так как T(k) = R(k)U(k) = U(k)**T( U(k)R(k) )U(k) = U(k)**T*T(k-1)*U(k) , то по индукции T(k) = (U(0)U(1)......U(k))**H*A(U(0)U(1)......U(k)). Следовательно, каждая матрица T(k) унитарно подобна матрице A. Не так заметно, но именно это центральная тема данного раздела, что матрицы T(k) почти всегда сходятя к верхней треугольной форме. То есть (7.3.2) почти всегда "сходится" к разложению Шура матриц A. Процесс (7.3.1) называют QR-итерациями, и он составляет основу наиболее эффективного алгоритма вычисления разложения Шура. Для того чтобы обосновать этот метода и установить характер его сходимости, рассмотрим сначала два других ирерационных процесса, представляющих самостоятельный интерес. 7.3.1. Степенной метод 7.3.2. Ортогональные итерации 7.3.3. QR-итерации 7.3.4. LR-итерации Приложение 7.4 Хессенбергова форма и вещественная форма Шура 324 В этом и следующем разделах мы покажем, как сделать QR-итерации (7.3.1) быстрым и эффективным методом вычисления разложения Шура. Поскольку большинство проблем собственных значений и собственных подпространств вещественные, мы остановимя на развитии вещественного аналога (7.3.1), который запишем следующим образом: H(0) = U(0)**H*A*U(0) for k = 1,2,.... H(k-1) = U(k)R(k) (QR-разложение) H(k) = R(k)U(k) end Здесь A -> R(n,n) , каждая матрица U(k) -> R(n,n) - ортогональная и каждая матрица R(k) -> R(n,n) - верхняя треугольнам. Трудность, связанная с этими вещественными итерациями, состоит в том, что матрицы H(k) могут никогда не "сойтись к строго "открывающей собственные значения" треугольный форме в случае, если матрица имеет комплексные собственные значения. По этой причине мы должны оставить бесплодные надежды и быть довольны вычислением альтернативного разложения, известного как вещественное разложения Шура. Заметим, что в численной линейной алгебре обычно, не нарушая общности, фокусируют внимание на вещественных матричных задач, поскольку большинство вещественных алгоритмов имеют очевидные комплексные аналоги. QR-итерации не являются исключением. Для того чтобы эффективно вычислять вещественное разложение Шура, мы должны тщательно выбрать начальное ортогональное преобразование подобия U(0) в (7.4.1). В чтастности, если мы выберем U(0) так, что матрица H(0) будет верхней хессенберговой, то количество работы на итерации уменьшается с O(n**3) до O(n**2). Начальное преобразование к хессенберговой форме (U(0)-вычисление) само по себе является очень важным алгоритмом и может быть реализовано последовательностью матричных операций Хаусхолдера. 7.4.1. Вещественное разложение Шура 7.4.2. Хессенбергов QR-шаг 7.4.3. Разложение Хессетберга 7.4.4. Трехуровневые варианты 7.4.5. Важные свойства матрицы Хессенберга 7.4.6. Сопровождающая матрица 7.4.7. Хессенбергово преобразование через преобразование Гаусса 7.5 Практический QR-алгоритм 334 Вернемся к хессенберговым QR-иреациям, которые мы запишем следующим образом: H(0) = U(0)**H*A*U(0) (Преобразование Хессенберга) for k = 1,2,.... H = UR (QR-разложение) H = RU end Цель этого раздела - описать сходимость матрицы H к верхней квазитреугольной форме и показать, как скорорсть сходимости можно увеличить, вводя сдвиги. 7.5.1. Исчерпывание 7.5.2. QR-итерации со сдвигами 7.5.3. Стратегия одинарного сдвига 7.5.4. Стратегия двойного сдвига 7.5.5. Стратегия неявного двойного сдвига 7.5.6. Полный процесс 7.5.7. Масштабирование 7.6 Методы вычисления инвариантных подпространств 343 Различные важные проблемы инвариантных подпространств могут быть решены, если найдено вещественное разложение Шура. В этом разделе мы обсуждаем как Вычислить собственные векторы, соответствующие некоторому подмножеству laymbda(A). Найти ортоганальный базис для заданного инвариантного подпространства. Блочно диагонализировать A, используя хорошо обусловленные преобразования подобия. Найти базис из собственных векторов, невзирая на их обусловленность. Найти приближенную каноническую форму Жордана матрицы A. Вычисление собственного вектора (инвариантного подпространства) разреженной матрицы обсуждается в другом месте. Смотри 7.3, а также гл. 8 и 9. 7.6.1. Выделение собственных векторов при помощи обратных итераций 7.6.2. Упорядочение собственных значений в вещественной форме Шура 7.6.3. Блочная диагонализация 7.6.4. Базис собственных векторов 7.6.5. Определение блочных жордановых структур 7.7 QZ-метод для Ax = lyambdaBx 353 Пусть A и B - две матрицы размера n x n. Множество всех матриц вида A - lambda*B c lambda -> C называют пучком. Собственные значения пучка - элементы множества lambda(A,B), определяемого как lambda(A,B) = {z -> C: det(A - zB) = 0 }. Если lambda -> lambda(A,B) и Ax = lambda*B*x, x != 0, (7.7.1) то вектор x называют собственным векотором пучка матриц (A - lambda*B). В этом разделе мы даем краткий обзор математических свойств обобщенной проблемы (7.7.1) и представляем устойчивый алгоритм ее решения. Важный случай, когда матрицы A и B - симметричные и матрицы B - положительно определенная, рассматривается в 8.7.2. 7.7.1. Основы теории 7.7.2. Обобщенное разложение Шура 7.7.3. Выводы о чувствительности 7.7.4. Хессенбергово-треугольная форма 7.7.5. Исчерпывание 7.7.6. QZ-шаг 7.7.7. QZ-процесс в целом 7.7.8. Вычисление обобщенных инвариантных подпространств Глава 8. Сииметрична проблема проблема собственных значений 368 Теория возмущения и алгоритмические разработки из предыдущей главы значительно упрощаются, когда матрица A симметричная. Симметричная пробема собственных значений с ее богатой математической структурой поистине является одной из наиболее притягательных проблем вычислительной алгебры. Мы начнем наше знакомство с короткого обсуждения математических результатов, лежащих в основе симметричной проблемы собственных значений. В 8.2 мы специализируем алгоритмы из 7.4 и 7.5 и получаем изящный симметричный QR-алгоритм. Вариант этой процедуры, позволяющей вычислять сингулярное разложение, подробно рассматривается в 8.3. Поскольку собственные значения симметричной матрицы могут быть определены различными способами, существует множесво QR-алгоритмов. Некоторые из этих методов описаны в 8.4. Одним из самых первых матричных алгоритмов, появившихся в литературе является метод Якоби для симметричной проблемы собственных значений. Интерес к этой процедуре возобновился благодаря параллельным вычислениям и поэтому мы отвели этому методу значительное место в 8.5. Другой высокопараллельный метод для проблемы собственных значений, вариант метода "разделяй и властвуй", рассмотрен в 8.6. Его сожно использовать для трехдиагональной проблемы, и он особенно интересен, так как успешно конкурирует с трехдиагонильными QR-итерациями даже в условиях однопроцессорных вычислений. В конце раздела мы рассматриваем A - lambda В проблему для важного случая, когда матрица A - симметричная, а матрица B - симмемтричная положительно определенная. Хотя не существует подходящих аналогов QZ-алгоритма для этой обобщенной проблемы собственных значений со специальной структурой, существует несколько удачных методов, которые могут быть использованы. Также обсуждается обобщенная проблема сингулярных значений. Большую часть материала этой главы можно найти в [Parlett, SEP]. 8.1 Математические основы 368 Симметричная проблема собственных значений имеет изящную и богатую теорию. Наиболее важные аспекты этой теориии рассматриваются в данном разделе. 8.1.1. Собственные значения симметричных матриц 8.1.2. Теорема о минимаксе и некоторые следствия 8.1.3. Дополнительные результаты 8.1.4. Чувствительность инвариантных подпространств 8.1.5. Закон инерции 8.2 Симметричный QR-алгоритм 376 Посмотрим теперь, как практический QR-алгоритм, полученный в гл. 7, может быть специализирован, если матрица A -> R(n,n) симметричная. Три замечания можно сделать сразу: Если матрица U(0)**T*A*U(0) = H верхняя хессенбергова, то H = T должна быть трехдиагональной. Сииметричность и трехдиагональная ленточная структуры сохраняются после выполнения QR- шага с одинарным сдвигом. Нет необходимости рассматривать комплексные сдвиги, так как lambda(A) => R. Эти урпощения вместе с приятными математическими свойствами симметричной проблемы собственных значений делаю алгоритмы этой главы очень привлекательными. 8.2.1. Преобразование к трехдиагональному виду 8.2.2. QR-итерации с явным одинарным сдвигом 8.2.3. Вариант с неявным сдвигом 383 192 8.3 Вычисление SVD 383 Существуют важные связи между сингулярным разложением матрицы A и разложениями Шура симметричных матриц A(**T)*A , A*A(**T) и [0 A(**T)] . В самом деле, если [A 0 ] U(**T)*A*V = diag(sigma1,....,sigma(n)) есть SVD матрицы A -> R(m,n) (m >= n), то V(**T)*(A(**T)*A)*V = diag(sigma(**2)(1),....,sigma(**2)(n)) -> R(n,n) (8.3.1) и U(**T)*(A*A(**T))*U = diag(sigma(**2)(1),....,sigma(**2)(n),0,....,0) -> R(m,m) (8.3.2) Кроме того, если U = [U(1)U(2)] n m-n и мы определим ортогональную матрицу Q -> R((m+n)x(m+n)) как [V V 0 ] Q = 1/root(2)[ ], [U(1) -U(1) root(2)U(2)] то [ 0 A(**T) ] Q(**T)[ ]Q = diag(sigma1,....,sigma(n),-sigma1,...,-sigma(n),0,....,0). [ A 0 ] m - n Эти связи с симметричной проблемой собственных значений возводяют нам приспособить математические и алгоритмические разработки из двух предыдущих разделов к проблеме сингулярных значений. 8.3.1. Теория возмущения и свойства 8.3.2. SVD-алгоритм 8.4 Некоторые специальные методы 392 Используя богатую математическую структуру симметричной проблемы собственных значений, можно придумать полезные альтернативы симметричному QR-алгоритму. Многие из этих методов предпочтительнее, когда требуется найти только несколько собственных значений и/или собственных векторов. Три таких метода описаны в этом разделе: бисекция, итерации с отношением Рэлея, ортогональные итерации с ускорением Ритца. 8.4.1. Бисекция 8.4.2. Итерация с отношением Рэлея 8.4.3. Ортогональные итерации с ускорением Ритца 8.5 Методы Якоби 399 Методы Якоби для сииметричной проблемы собственных значений в настоящее время привлекает к себе внимание, потому что они являются существенно параллельными. Во время их работы выполняется последовательность ортогонально подобных обновлений A <- Q(**T)*A*Q, обладающих тем свойством, что каждая новая матрица A, даже плотная, является "более диагональной", чем ее предшественник. В конце концов, внедиагональные элементы становятся достаточно малыми для того, чтобы объявить их равными нулю. После обзора основных идей, лежащих в основе метода Якоби, мы разработаем параллельную процедуру Якоби, которая годится для кольцевой многопроцессорной системы. Все процедуры Якоби для симметричной проблемы собственных занчений из этого раздела имеют SVD-аналоги. 8.5.1. Идея метода Якоби 8.5.2. Симметричное разложение Шура размера 2 x 2 8.5.3. Обновления в методе Якоби 8.5.4. Классический алгоритм Якоби 8.5.5. Алгоритм циклический по строкам 8.5.6. Барьерный метод Якоби 8.5.7. Анализ ошибок округления 8.5.8. Сравнение с симметричным QR-алгоритмом 8.5.9. Параллельное упорядочение 8.5.10. Кольцевая процедура 8.5.11. Блочная процедура Якоби 8.5.12. SVD-процедура Якоби 207 8.6 Метод разделяй и властвуй 412 8.6.1. Расщепление 8.6.2. Объединение разложений Шура 8.6.3. Собственная система матрицы D + pzz(t) 8.6.4. Практический синтез 8.6.5. Полный процесс с параллелизмом 210 8.7 Более общие проблемы собственных значений 418 8.7.1. Математические основы 8.7.2. Методы для симметрично-определенной проблемы 8.7.3. Обобщенная проблема сингулярных значений Глава 9. Методы Ланцоша 426 В этой главе мы изучаем метод Ланцоша - метод решения некоторого класса больших разреженных симметричных спектральных задач Ax = lambda x. В методе используются частичные трехдиагонализации заданныой матрицы A. Однако в отличие от подхода Хаусхолдера здесь не строятся какие-либо промежуточные полные матрицы. В равной степени важно то, что информация об экстремальных собственных значениях для A обычно появляется задолго до построения полной трехдиагонализации. Это делает методв Ланцоша особенно полезным в тем случаях, когда для A требуется найти небольшое число наибольших или наименьших собственных значений. В 9.1 представлены вывод и свойства метода в точной арифметике. Детально рассмотрены ключевые аспекты теории Каниэля-Пейджа. Это теория объясняет необычные свойства сходимости процесса Ланцоша. К несчастью, практическое использование метода Ланцоша несколько затрудняется ошибками округления. Центральная проблема - это потеря ортогональности получаемых итерационно векторов Ланцоша. С этим можно справиться нескольким способами - они обсуждаются в 9.2. В заключительном разделе мы показываем, как "идея Ланцоша" может быть применена к различным задачам, связанным с сингулярными числами, наименьшими квадратами и линейными уравнениями. Обсуждаются также и несимметричный процесс Ланцоша. В 9.3 особый интерес предстваляет разработка метода сопряженных градиентов для симметричных положительно определенных линейных систем. Изучение связи метода Ланцоша и метода сопряженных градиентов будет продолженио в следующей главе. 214 9.1 Выводы свойства сходимости 426 Пусть имеется большая разреженная симметричная матрица A -> R(n,n) и нас интересует небольшое число ее наибольших и (или) наименьших собственных значений. Эта задача может быть решена с помощью метода, приписываемого Ланцошу(1950). Метод строит последовательность трехдиагональных матриц T(j) с тем свойством, что экстремальные собственные значения для T(j) -> R(j,j) с ростом j дают все более точные оценки экстремальных собственных значений для A. В этом разделе мы выводим метод и исследуем его свойства в точной арифметике. 9.1.1. Подпространства Крылова 9.1.2. Трехдиагонализация 9.1.3. Окончание и оценки погрешностей 9.1.4. Теория сходимости Каниэля-Пейджа 9.1.5. Сравнение степенного метода с методом Ланцоша 9.1.6. Сходимость внутренних собственных значений 217 9.2 Практические процедуры Ланцоша 433 На поведение итераций Ланцоша сильнейшее воздействие оказывают ошибки округления. Основная трудность вызвана потерей ортогональности среди векторов Ланцоща. Это явление затуманивает вопрос об окончании процесса и усложняет соотношения междц собственными значениями матрицы и трехдиагональных матриц . Трудности, навеянные этим явлением, вместе с созданием совершенно устойчисовго метода хаусхолдеровой трехдиагонализации, объясняют, почему численные аналитики не обращали внимания на алгоритм Ланцоша в течении 50-х и 60-х годов. Однако интерес к методу оживили как теория Каниэля-Пейджа, так и необходимость решать большеие разреженные спектральные задачи, размер которых увеличивался вместе с производительностью компьютеров. Поскольку для получения хороших приближений к экстремальныхм собственным значениям обычно требуется намного меньше n итераций, метод Ланцоша привлекает скорее не как соперник подхода Хаусхолдера, а как средство для работы с разреженными матрицами. Успешные реализации метода Ланцоша включает в себя гораздо больше, чем только закодированные соотношения (9.1.3). В этом разделе мы дадим представление о некоторых практических идеях, предложенных для того, чтобы сделать процедуру Ланцоша работоспособной. 9.2.1. Реализация в точной арифметике 9.2.2. Ошибки округления 9.2.3. Метод Ланцоша с полной переортогонализацией 9.2.4. Выборочная ортогонализация 9.2.5. Проблема теневых собственных значений 9.2.6. Блочный Ланцош 9.2.7. s-Шаговый алгоритм Ланцоша 222 9.3 Приложения и обобщения 442 В этом разделе мы вкратце расскажем о том, как итерации Ланцоша могут быть приспособлены для того, чтобы решать большие разреженные линейные системы и задачи наименьших квадратов. Мы обсудим также процессы Арнольди и несимметричного метода Ланцоша. 9.3.1. Симметричные положительно определенные системы 9.3.2. Симметричные неопределенные системы 9.3.3. Двухдиагонализация и SVD 9.3.4. Наименьшие квадраты 9.3.5. Идея Арнольди 9.3.6. Несимметричная трехдиагонализация Ланцоша Глава 10. Итерационные методы для линейных систем 452 В конце предыдущей главы было показано, каким образом итерации Ланцоша можно использовать для решения различных линейных уравнений и задач наименьших квадратов. Разработанные методы подходят для больших разреженных задач, так как они не требуют разложения соответствующей матрицы. В этой главе мы продолжим обсуждение процедур решения линейных уравнений, обладающих этим же свойством. Первый раздел - это беглый обзовр классических итерационных методов: Якоби, Гаусса-Зейделя, SOR, чебышевских полуитераций и т.д. Мы даем краткое изложение этих методов, потому что наша главная цель в этой главе - основательно представить метод сопряженных градиентов. В 10.2 мы подробно излагаем эту важную технику, естественным образом отталкиваясь от метода наисеорейшего спуска. Напомним, что в 9.3 метод сопряженных градиетов уже появился в связи с итерациями Ланцоша. Метод выводится заново для того, чтобы объяснить некоторые его практические варианты, которым посвящен 10.3. Мы предупреждаем читателя о непоследовательности в обозначениях этой главы. В 10.1 разработка методов на "i,j-уровне" делает необходимым использование верхних индексов: x(k)(i) обозначает i-ю компоненту вектора x(k). В других разделах, однако, развитие алгоритмов происходит без явного упоминания компонент векторов или матриц. Поэтому в 10.2 и 10.3 мы можем отказаться от верхних индексов, а последовательности векторов обозначать как {x(k)}. 227 10.1 Стандартные итерации 452 В гл. 3 и 4 процедуры решения линейных уравнений были связаны с разложением матрицы коэффициентов . Методы такого типа называются прямыми методами. Если матрица A большая и разреженная, то искомые составляющие разложений могут быть плотными, и это ограничивает практическое применение прямых методов. Одним из исключений является случай ленточной матрицы A (см. 4.3). Но даже во многих задачах с летночными матрицами сама лента оказывается разреженной, внося трудности в реализацию алгоритмов типа ленточного метода Холецкого. Одна из причин огромного интереса к процедурам решения разреженных линейных уравнений связана с тем значением, которое придается умению находить численные решения уравнений в частных производных. В действительности именно исследователи в области вычислительных методов для уравнений в частных производных ответственны за многие из подходов к разреженным матрицам, которые теперь используются повсеместно. Грубо говоря, есть два подхода к разреженной задаче Ax = b. Один из них заключается в выборе подходящего прямого метода и его адаптации с тем, чтобы учесть разреженность в A. Типичные стратегии адаптации связаны с разумным использованием структур данных и со специальными, минимизирующими заполнения стратегиями выбора ведущего элемента. Есть обширная литература по этим вопросам; заинтересованному читателю следует обратиться к книгам George, Liu (1981) и Duff, Erisman, Reid (1986). Противоположностью прямым методам являются итерационные методы. Эти методы порождают последовательность приближенных решений {x(k)}, и по существу A участвует лишь в матрично-векторных умножениях. При оценивании качества итерационного метода неизменно в ценре внимания вопрос о том, как быстро сходятся итерации x(**k). В этом параграфе мы представим некоторые основные итерационные метода, обсудим их практические реализации и докажем несколько характерных теорем о поведении итераций. 10.1.1. Итерации Якоби и Гаусса-Зейделя 10.1.2. Расщепление и сходимость 10.1.3. Практическая реализация Гаусса-Зейделя 10.1.4. Последовательная верхняя релаксация 10.1.5. Метод Чебышевских полуитераций 10.1.6. Симметричный метод SOR 231 10.2 Методы сопряженных градиентов 461 Трудность, связанная с SOR, чебышевскими полуитерациями и такого же типа методами, заключается в том, что они зависят от параметров, правильный выбор которых иногда бывает затруднителен. Например, для того чтобы чебышевское ускорение было успешниым, нам нужны хорошие оценки для наибольшего и наименьшего собственных значений соответствующей итерационной матрицы M(**-1)N. Если эта матрица не устроена по-особому, то полуыение их в аналитическом виде, скорее всего, невозможно, а вычисление дорого. В этом разделе мы представим метод, у котором нет этой трудности,- хорошо известный метод сопряженных градиентов Хестенса-Штифеля. Мы вывели этот метод в 9.3.1 при обсуждении - как приспособить алгритм Ланцоша к решению симметричных положительно определенных линейных систем. Здесь мы выводим метод сопряженных градиентов заново для того, чтобы подготовить почву для 10.3, где будет показана его связь с другими итерационными схемами для Ax = b и будут описаны некоторые из полезных его обобщений. 10.2.1. Наискорейший спуск 10.2.2. Произвольные направления спуска 10.2.3. A-сопряженные направления спуска 10.2.4. Метод сопряженных градиентов 10.2.5. Несколько крайне необходимых наблюдений 10.2.6. Связь с алгоритмом Ланцоша 10.2.7. Некоторые практические детали 10.2.8. Сходимость метода 236 10.3 Сопряженные градиенты с предобусловливанием 471 Предыдущий параграф мы закончили замечанием о том, что метод сопряженных градиентов работает хорошо для матриц, которые либо хорошо обусловлены, либо имеют лишь небольшое число различных собственных значений. (Последнее относится к случаю, когда A есть малоранговое возмущение единичной матрицы). В этом параграфе мы покажем, как проводится предобусловливание линейной системк с тем, чтобы матрицы коэффициентов приобретала одну из этих прекрасных особенностей. Наше изложение будет весьма сжатым и неформальным. Более основательно эти вопросы представлены в Golub, Meurant (1983) и Axelsson (1985). 10.3.1. Вывод 10.3.2. Предобусловливатели, связанные с неполным разложением Холецского 10.3.3. Неполные блочные предобусловливатели 10.3.4. Идеи декомпозиции области 10.3.5. Полиноминальные предобусловливатели 10.3.6. Заключительное предложение Глава 11. Функции от матриц 482 Задача вычисления функции f(A) от матрицы A размера n x n часто встречается в теории управления и других прикладных областях. Грубо говоря, если скалярная функция f(z) определена на спектре lambda(A), то функцию от матрицы f(A)определяют, подставляя матрицы A вместо переменноой z в "формулу" для f(z). Например, если f(z) = (1 + z)/(1 - z) и 1 не принадлежит lambda(A), то f(A) = (I + A)/(I - A). Эти вычисления становятся особенно интресными, когда функция f трансцендентная. Один из подходов в этой более сложной ситуации состоит в вычислении какого-нибудь спектрального разложения A = YB/Y и использовании формулы f(A) = Yf(B)/Y. Если матрица B достаточно простая, то часто можно вычислить f(B) непосредственно. Это показано в 11.1 для разложения Жордана и Шура. Не удивительно, что на основе последнего разложения получается более устойчивая f(A) процедура. Другой подход к вычислению функции от матрицы состоит в аппроксимации требуемой функции f(A), леко вычисляемой функцией g(A). Например, функция g может быть усеченным рядом Тейлора, аппроксимирующим функцию f. Оценки погрешности, связанной с такой аппроксимацией функции от матрицы, даны в 11.2. В последнем разделе мы обсуждаем специальную и очень важную задачу вычисления матричной экспоненты e(A). 242 11.1 Спектральные методы 482 Для заданной матрицы A размера n x n и скалярной функции f(z) существует несколько способов определить функции от матрицы f(A). Неформальный способ определения - подставить матрицу вместо переменной z в формулу для f(z). Например, если p(z) = 1 + z и r(z) = (1+z/2)/(1-z/2) при z != 2, то, конечно, разумно определить p(A) и r(A) как p(A) = I + A и r(A) = (I + A/2)/(I-A/2), 2 !-> laymbda(A). Подстановка A вместо z также срабатывает и в случае трансцендентных функций, например, e(**A) = SUM(k=0,00)(A**k)/k!. Однако для того чтобы сделать последующие алгоритмические разработки корректными, мы должны более аккуратно определить f(A). 11.1.1. Определение Существует много способов строгого определения понятия функуции от матрицы. Смотри [Rinehart, 1955]. Возможно, наиболее изящный подход основан на применении контурного интеграла. Пусть функция аналитична на замкнутом множестве, ограниченном замкнутым контуром , который окружает спектр . Определим как матрицу В этом определении можно сразу узнать матричный вариант интегральной теоремы Коши. Интеграл определяется поэлементно Заметим, что элементы матрицы аналитичны на и матрица опредене в тех случаях, когда функция аналитична в окрестности . 11.1.2. Жорданова характеризация Определение , совершенно бесполезное с вычислительной точки зрения, можно использовать для вывода более практичных характеризаций . Например, если матрица определена, и 11.1.3. Метод разложения Шура 11.1.4. Метод блочного разложения Шура 245 11.2 Аппроксимационные методы 488 Теперь мы рассмотрим класс методов вычисления функций от матрицы, в которых собственные значения не играют главную роль. Эти методы основаны на следующей идее: если функция g(z) аппроксимирует f(z) на спектре laymbda(A), то матрица g(A) аппроксимирует матрицу f(A), например: e(**A) = I + A + A**2/2! + .... + A**q/q! . Мы начнем с оценки нормы ||f(A) - g(A)|| при помощи представлений Жордана и Шура функции от матрицы. И продолжим некоторыми замечаниями о величине матричных многочленов. 11.2.1. Жорданов анализ 11.2.2. Анализ на основе разложения Шура 11.2.3. Тейлоровы апроксиманты 11.2.4. Оценка многочлена от матрицы 11.2.5. Вычисление степеней матрицы 11.2.6. Интегральные функции от матрицы 248 11.3 Матричная экспонента 495 Одной из наиболее часто вычисляемых функций от матрицы является экспонента e**At = SUM(k=0,00)(At**k)/k!. Было предложено множество алгоритмов вычисления матрицы , но, как указано в обзорной статье Молера и Ван Лоана [Moler, Van Loan, 1978], большинство из них обладает сомнительными вычислительными свойствами. Для того чтобы пояснить, какие имеются вычислительные трудности, мы проведем краткий анализ возмущения матричной экспоненты и затем используем его для оценки одного из лучших e**At алгоритмов. 11.3.1. Теория возмущений 11.3.2. Метод аппроксимации Паде 11.3.3. Некоторые выводы об устойчивости Глава 12. Специальные разделы 501 В последней главе мы обсуждаем набор задач, в котором представлены важные приложения метода SVD и QR-разложения. Сначала мы рассмотрим минимизацию в смысле наименьших квадратов с ограничениями. В 12.1 рассматриваются два типа ограничений - квадратичные неравенства и линейные равенства. Следующие два раздела связаны с вариантами стандартной LS-задачи. В 12.2 мы рассмотрим, как может быть аппроксимирован вектор b некоторым подмножеством столбцов матрицы A. Этот способ иногда применяют, если матрицы A неполного ранга. В 12.3 мы рассматриваем вариант обычной регрессии, известной как общая задача наименьших квадратов, которая возникает, если матрица A искажена ошибками. Другие приложения метода SVD обсуждаются в 12.4, где рассмотрены различные подобласти вычислений. В 12.5 обсуждаются некоторые варианты симметричной проблемы собственных значений. В 12.6 мы исследуем изменения QR-разложения, когда к матрице A = QR добавляется матрица ранга один. 251 12.1 Задача наименьших квадратов с ограничениями 501 Задача наименьших квадратов иногда естественно рассматривать как минимизацию величины ||Ax - b||(2) на некотором подмножестве R(n). Например, мы хотим определсть вектор b, как наилучшее приближение для Ax, с ограничением, что x является единичным вектором. Или, возможно, решение определяет подходящую функцию f(t), которая описывается значенияи в конечном числе узлов. Это может приветси к задаче наименьших квадратов с ограничениями типа равенств. В этом разделе мы покажем, как эти задачи могут быть решены с использованием QR-разложения и SVD. 12.1.1. Задача LSQI 12.1.2. LS-минимизация на сфере 12.1.3. Гребневая регрессия 12.1.4. Задача наименьших квадратов с ограничениями типа равенств 12.1.5. Метод взвешивания 255 12.2 Выбор подмножеств при помощи SVD 509 Как упоминалось в 5.5, LS-задча неполного ранга min||Ax - b||(2) может быть заменена аппроксимацией решния минимальной нормы x(ls) = SUM(i=1,r)( (u(i)(**T)*b)*v(i)/sigma(i) ), r = rank(A), с вектором x(r~) = SUM(i=1,r~)(u(i)(**T)*b(i))/sigma(i), r~ <= r, где A = U*Sig*V(**T) = SUM(i=1,r)(sigma(i)*u(i)*v(i)**T). SVD-разложение матрицы A и (r~)-вычисленная оценка для r. Отметим, что вектор x(r~) минимизирует ||A(r~)*x -b||(2) , где матрица A(r~) = SUM(i=1,r)(sigma(i)*u(i)*v(i)**T) является приближением матрицы A ранга r~. См. теорему 2.5.2. Замена матрицы A на матрицу A(r~) в LS-задаче выполняется для фильтрации малых сингулярных значений и имеет большой смысл в ситуациях, когда матрица A получаетя из возмущенных данных. Однако в других приложениях неполнота ранга означает избыточность факторов, которые содержит лежащая в основе модель. В этом случае для создателя модели может быть неинтересным такое приближение, как A(r~)x(r~), которое учитывает все n избыточных факторов. Вместо этого можно найти приближение Ay, где y имеют около r~ ненулевых компонент. Позиции ненулевых элементов вектора определяют некоторые стобцы матрицы A, т.е. некоторые факторы модели, которые используются для приближения вектора наблюдений b. Задача выбора этих столбцов называется выбором подмножества и составляет содержание данного раздела. 12.2.1. Метод QR со столбцовым выбором 12.2.2. Использование метода SVD 12.2.3. Еще раз о противоречии между независимостью столбцов и невязкой 258 12.3 Общая задача наименьших квадратов 514 Задача миниминзации ||D(Ax-b)||(2), где матрица A -> R(m x n), а матрица D = diag(d1,...dm) является невырожденной, может быть переформулирована следующим образом: min ||Dr||(2), r -> R(m) b+r -> range(A) В этой задаче по умолчанию существует предположение, что ошибки заключены в векторе "наблюдений" b. Когда ошибка также присутствует в "данных" матрицы A, то более естественно рассмотреть задачу min ||D[Er]*T||(F)*E -> R(mxn), r -> R(m), (12.3.2) где матрицы D = diag(d1,...dm) и T = diag(t1,...,t(n+1)) невырождены. Эта задача, рассмотренная в работе Golub и Van Loan, была названа общей задачей наименьших квадратов (TLS). Если минимизатор [E(0)r(0)] может быть найден для задачи (12.3.2), то любой вектор x, удовлетворяющий системе (A + E(0))*x = b + r(0), называется TLS-решением. Однако следует понимать, что задача (12.3.2) может не иметь совместного решения. Например, если [1 0] [1] [0 0] A = [0 0], b =[1], D = I(3), T = I(3), E(e)=[0 e], [0 0] [1] [0 e] то для всех e > 0 b -> range(A + E(e)). Тем не менее не существует наименьшего значения ||[E*r]||(F), для которого b + r -> range(A + E). Обобщение задачи (12.3.2) получается, если мы допускаем наличие многих правых частей. В частности, если матрица B -> R(mxk), то мы имеем задачу min ||D[ER]T||(F). range(B+R) -> range(A+E) где E -> R(mxn) и R -> R(mxk), а матрицы D = diag(d1,....,dm) и T = diag(t1,...,t(n+k)) неыврождены. Если [E(0)R(0)] разрешает задачу (12.3.3), то любой X -> R(mxk), удовлетворяющий системе (A + E(0))X = (B + R(0)), называется TLS-решением для задачи (12.3.3). В этом разделе мы рассмотрим некоторые из математических свойств общей задачи наименьших квадратов и покажем, как она может быть решена с использованием SVD. 12.3.1. Математическая основа 12.3.2. Вычисления для случая k = 1 12.3.3. Геометрическая интерпретация 260 12.4 Сравнение подпространств при помощи SVD 518 Иногда необходимо исследовать связь между двумя заданными подпространствами. Насколько они близки? Пересекаются они или нет? Может ли одно быть "повернутым" в другое? И так далее. В этом разделе мы покажем, как можно ответить на подобные вопросы при помощи сингулярного разложения. 12.4.1. Поворот подпространств 12.4.2. Пересечение ядер 12.4.3. Углы между подпространствами 12.4.4. Пересечение подпространств 262 12.5 Модифицированные задачи на собственный значения 523 В этом разделе рассматриваются некоторые варианты стандартной задачи на собственные значения. В обсуждении широко используются матричные методы, описанные в этой книге. 12.5.1. Стационарные значения квадратичной формы с ограничениями 12.5.2. Обратная задача на собственные значения 12.5.3. Одноранговая модификация задачи на собственые значения 12.5.4. Вторая обратная задача на собственные значения 265 12.6 Модификация QR-разложения 528 Во многих приложениях требуется заново вычислить разложение заданной матрицы A -> R(mxn) после того, как матрица была незначительно изменена. Например, нам задано QR-разложение матрицы A и требуется вычислить QR-разложение матрицы A, которая получается (а) прибавлением произвольной одноранговой матрицы к A, (б) добавлением строки (столбца) к матрице A, (в) удалением строки (столбца) из матрицы A. В данном разделе мы покажем, что в ситуациях, подобных этой, гораздо эффективнее "модифицировать" QR-разложение матрицы A, чем выполнить его заново. Перед началом упомянем, что существуют также способы для модификаций разложений PA = LU, A = GG(**t) и A = LDL(**T). Однако эти разложения могут быть довольно чувствительны к модификациям, потому что требуется выбор ведущего элемента и, кроме того, когда мы изменяем положительно определенную матрицу, то результат может не быть положительно определенным. См. Gill, Golub, Murray и Saunders (1974) и Stewart(1979c). Следуя этим результатам, мы кратко обсудим гиперболические преобразования и их использование в задаче модификации разложения Холецкого после удаления строки. 12.6.1. Одноранговые модификации 12.6.2. Добавление или удаление столбца 12.6.3. Добавление или исключение строки 12.6.4. Методы гиперболического преобразования Предметный указатель Алгебраическая кратность (algebraic multiplicity) 305 Алгоритмы: (algorithms) Аазена (Aassen's) 155 барьерный метод Якоби (threshold Jacobi) 403 бисекция (bisection) 393 блочная циклическая редукция (block cyclic reduction) 162-164 блочный Ланцоша (block Lanczos) 438 быстрое вращение Гивенса (fast Givens rotation) 210 возведение в квадрат exp(A) (squaring) 498, 499 Гаусса-Зейделя (Gauss-Seidel) 454 двоичное возведение в степень (binary powering) 493 двухдиагонализация Хаусхолдера (Hausholder bidiagonalization) 218 диагонального выбора (diagonal pivoting) 156 Дулитла (Doolittle reduction) 99 жорданово разложение (Jordan decomposition) 350 исключение Гаусса (Gauss elimination) блочный (block) 113 полный выбор ведущего элемента (complete pivoting) 114 с внешним произведением (outer product) 96 с частичным выбором: версия с внешним произведением (outer product with pivoting) 109 с частичным выбором: gaxpy-версия (gaxpy with pivoting) 111 gaxpy-версия LU-разложения (gaxpy LU) 97 итерация с разложением Релея (Rayleigh quotient eteration) 440 итерация Якоби (Jacobi iteration) 453 классический метод Грамма-Шмидта (classical Gram-Schmidt) 201 ланцошева двухдиагонализация (Lanczos bidiagonalization) 446, 447 Левинсона (Levinson) 175 ленточный метод Холецкого (band Cholesky) 146 леточная обратная подстановка (band backward substution) 143 ленточная прямая подстановка (band forward substitution) 143 матричная операция gaxpy (matrix-matrix gapxy operation) 22 метод Арнольди (Arnoldi method) 448 метод Бартельса-Стюарта (Bartels-Stewart nethod) 486 метод оценки обусловленности (condition estimator) 124 метод gapxpy с обратным ходом (backward gaxpy) 458 неявный симметричный QR-шаг со сдвигом Уилкинсона (implicit symmetric QR step with Wilkinson) 380 нормальные уравнения (normal equations) 207 обратная подстановка (backward substitution) 88, 89 ортогонализация Грама-Шмидта (Gram-Schmidt orthogonalization) 202 ортогональные итерации (orthogonal iterations) 318 Партлета-Рида (Partlett-Reid) 151 параллельная кольцевая факторизация (parallel ring factorization) 277 пересечение подпространств (intersection of subspaces) 522 пересечение ядер (null-space intersection) 519 последовательная верхняя релаксация (sequential over-relaxation SOR) 456 процедура Ланцоша (Lanczos) 433 прямая подстановка (forward substitution) 87, 89 прямое накопление (Hausholder matrix accumulation) 184 решение блочно-диагональной системы (block diagonal system solving) 160-161 решение задачи Прокруста (Procrust) 518 решение симметричной трехдиагональной положительно определенной системы (positive definite tridiagonal system solver) 147 решение системы Вандермонда (Vandermond system solving) 169 симметричный метод (symmetric successive over-relaxation) 459 скалярное произведение (dot product) 18 скалярное произведение матриц (matrix-matrix dot) 23,25 сопряженные градиенты с предобусловливанием (preconditioned conjugate gradients) 472,473 стационарные значения квадратичной формы с ограничениями (starionary values with constraints) 523 степенные итерации (power iterations) 316 Тренча (Trench) 177 углы между подпространствами (angles between subspaces) 521 умножение матриц с использованием внешних произведений (matrix-matrix outer product) 25,26 умножение матрицы на вектор (matrix-vector row) 19 умножение на матрицу Гивенса (Givens rotation times matrix) 188 умножение на матрицу Хаусхолдера слева (Hausholder reflection times matrix) 184 умножение треугольных матриц (triangular matrix multiplication) 30 функция от треугольной матрицы (function of triangular matrix) 486 Хаусхолдера трехдиагонализация (Hausholder triangularization) 377 Хаусхолдера QR-разложение с выбором верхнего столбца (Haushilder QR) 216 хесенбергово-треугольное преобразование (Hessenberg-trangular reduction) 356-358 Холецкого параллельная реализация (Cholesky parallel) 273,275,284,285 Холецкого разложение А-хВ (Cholesky decomposition) 420 Холецкого сопряженных направлений (conjugate gradients) 445,467 циклический метод Якоби (cycllic Jacobi) 402 чебышевские полуитерации (Chebyshev semi-iterative) 457 Штрассена (Strassen) 43 экспонента матрицы (matrix exponentiation) 559 gaxpy 22 LDLt- разложение () 128,130,131 LDMt- разложение () 127,129 LU- разложение матрицы Хессенберга (LU-Hessenberg) 145 QR-несимметричный (QR unsymmetric) 340 QR-симметричный (QR symmetric) 380,381 QR-итерации с одинарным сдвигом (single shift QR iteration) 336 QR-разложение: преобразование Гивенса (QR Givens) 198 QR-разложение: преобразование Хаусхолдера (QR Givens) 196 QR-шаг Френсиса (Francis QR step) 340 QR-шаг (QR step) 361-362 QZ-шаг (QZ step) 440 S-шаг (S-step Lanczos) 440 saxpy 18 SVD 390 SVD-процедура Якоби (SVD procedure Jacobi) 408,409 SVD-шаг Голуба-Кахана (SVD Golub-Kahan step) 389 TLS 516 Анализ ошибок (error analysis) -- обратный (backward) 70 -- прямой (forward) 70 Аппроксимация матричной функции (approximation of a matrix function) 488 Арифметика с округлением (rounded arithmetic) 66 - с отбрасыванием разрядов (shopped arithmetic) 66 База (base) 66 Базис (basis) 56 Бисекция (bisection) 393 Блочная матрица (block matrix) 36,37 Блочное диагональное доминирование (block diagonal dominance) 161 Блочные алгоритмы (block algorithms) 36,53 двухдиагональные системы (bidiagonal systems) 161 Ланцоша (Lanczos) 438 повторное использование данных (data re-use) 53 тридиагональные системы (tridiagonal systems) 159 циклическая редукция (cyclic reduction) 163,164 Якоби (Jacobi) 407,408 Быстрое преобразование Гивенса (fast Givens transformation) 192,193 Вектор (vector) - Гаусса (Gauss) 93 - Хаусхолдера (Hausholder) 183 Векторно-конвейерные вычисления (vector pipeline computing) 45 Векторно-конвейерный компьютер (vector pipeline computer) 45 270 Векторные вычисления (vector computing) 46,47 -- конвейерные (pipelining) 40,46,47 -- регистры (register) 47 - нормы (vector norms) 59,60 - обмены (vector touch) 52 - обозночения (notations) 18 - операции (operations) 18 Векторы Ланцоша 428 - Шура (Schur vectors) 303 Взаимные перестановки (interchange permutations) 107 Взвешивание по столбцам (column weighting) 230 -- строкам (row weighting) 230 Внедиагональный элемент (off-diagonal elements) 399 Внешниее произведение (outer product) 17 -- блочное (block) 40 Вращение Гивенса (Givens rotations) 187,188,198,199 - Якоби (Jacobi) 400 Выбор ведущего элемента (pivoting) в методе Аазена (Aasen) 155 --- полный (complete) 114 --- симметричный (symmetric matrices) 139,157 ---, столбец (column) 214,215 --- частичный (partial) 108 - подмножества (subset selection) 509 Выборочная ортогонализация (selective orthogonalization) 437 Вычисление ортонормированного базиса (orthonormal basis computation) 203 - скалярного произведения с накоплением (dot product accumulation) 69 --- погрещности округления (roundoff errors) 66 Гауссово исключение (Gaussian elimination) 23 -- ошибки округления (roundoff errors) 102 Геометрическая кратность (geometric multiplicity) 305 Гиперболическое преобразование (hyperbolic transformations) 532,533 Главные векторы и углы между подпространствами (principal angles and vectors) 521 Гребневая регрессия (ridge regression) 504 Двоеточие (colon) 21,32 Двоичное возведение в степень (binary powering) 493 Двойная точность (double precision) 69 Двухдиагонализация (bidiagonalization) 218 - метод Ланцоша (Lanczos) 446 - с приведением матрицы к верхнедианональному виду (upper triangularizing first) 219 Двухдианональная матрица (bidiagonal matrix) 30 Декомпозиция области (domain decomposition) 475, 476 Детерминант (determinant) 51,52,58 -, вырожденность (singularity) 83 Диагональная форма (diagonal form) 305 Динамическое распределение работы (dynamically sheduled algorithms) 257,261 Динамический резерв задач (pool of task) 256,270 Дифференцирование матрицы (differentiation) 58 Длина вектора (vector length) 47 Добавление или исключение строки (row addition or deletion) 531 Доминирующее собственное значение (dominant eigenvalue) 317 Доминирующий собственный вектор (dominant eigenvector) 317,318 Дополнение Шура (Schur comlement) 101 Дробление вычислений (granularity) 247,264,265 -- систолическая модель (systolic model) 266 Единичная матрица (identity matrix) 57 - ошибка (unit roundoff) 66 Единичный шаг выборки (unit stride) 50 Жорданово представление функции от матрицы (Jordan decomposition of matrix function) 488 - разложение (decomposition) 305,306 -- вычисление (computation) 350 Жордановы блока (Jordan block) 306 Задача наименьших квадратов (least square (LS) problem) 221 --- невязка (resudual) 206 --- неполного ранга (rank deficient) 221,222 --- основные решения (basic solutions) 224 --- полного ранга (full rank) 224 --- решение с минимальной нормой (minimum norm solution) 222 --- с ограничениями (constrained least squares) 501,503 ---- типа равенств (equality constrained LS) 505 --- чувствительность к возмкщениям (perturbation) 223 - Прокруста (Procrust problem) 518,519 -Юла-Уолкера (Yule-Walker problem) 172 Закон Амдаля (Amdahl's law) 52 - инерции Сильвестра (Sylvester law of inertia) 374 --- теорема (?) 374 Защищенные переменные (protected variables) 258 Иерархическая память (hierarchical memory) 49 Инвариантное подпростраство (invariant subspace) 300, 305 -- доминирование (dominating) 318 -- вектор Шура (Schur vector) 303 -- возмущения (perturbation) 371 -- прямая сумма (direct sum of) 305 -- чувствительность (sensitivity) 371, 372 Индекс профиля (profile index) 149 Инерция симметрической матрицы (enertia of symmetric matrix) 374 Инициализация (initialization) 258 Интервал значений показателя (exponent range) 66 Исчерпывание (deflation) 335 - хессенбергово-треугольной формы (Hessenberg-triangular form) 358 - QR-алгоритм (QR) 335 Итерации Гаусса-Зейделя (Gauss-Seidel iterations) 507 --- предобуславливатель (preconditioner) 477, 478 --- при решении уравнения Пуассона (solving Poisson equation) 456 - с одинарным шагом (single shift iterations) 377 -- отношением Релея (Rayleigh quotient iterations) 395 ---- QR-алгоритм (QR) 396 ---- симметрично-определенный пучок (symmetric-definite pencil) 421 Итерационная матрица (iteration matrix) 455 Итерационное уточнение (iterative improvement) 232 -- для метода наименьших квадратов (for LS) 232 -- линейных систем (for linear systems) 121 122 -- с фиксированной точностью (fixed precision) 123 -- со смешанной точностью (mixed precision) 122 Итерационные методы (iterative methods) 452 Итерация Якоби для метода (Jacobi eiteration fof SVD) 408 409 --- симметричной задачи собственных значений (symmetric eigenproblem) 399 Квадратичная форма (quadratic form) 523 Квадратный корень матрицы (squfre root of a matrix) 494 Классические итерации Якоби для собственных значений (classical Jacobi iterations for eigenvalues) 401 Кольцевые процедуры (ring algorithms) 273 -- другие (others) 276 -- Холецкого (Cholessky) 273, 274 -- Якоби (Jacobi) 405 Кольцо (ting) 240 Комплексная матрица (complex matrix) 27 Ковейеризация (pipelining) 45 - сложения (addition) 45, 46 - saxpy 46 Конфлюэнтные системы Вандермонда (confluent Vandermonde matrix) 170 Косинус матрицы (consine of a matrix) 490 Кратнасть алгебраическая (algebraic multiplicity ) 305 - геометрическая (geometric) 305 - собственного значения (multiplicity of eigenvalues) 305 Кратные собственные значения (multiple eigenvalues) 305 --- матричная функция (?) 486 --- триангуляция Ланцоша (?) 439 Критический разрез алгоритма (critical section) 257 Кросс-валидации (cross-validation) 504 Круговое распределение (wrap mapping) 274 Крылова матрица (Krylov matrix) 330 - подпространство (Krylov subspace) 426, 427 Ленточное LU-разложение (band LU factorization) 142, 143 Ленточные алгоритмы (band algorithms) исключение Гаусса (Gauss elimination) 144 ленточный метод Холецкого (band Cholesky) 146 решение треугольных систем (triangular system) 142 LU-разложение матрицы Хессенберга (LU-Hessenberg) 144, 145 Ленточный метод Холецкого (band Cholesky) 146 - параллельный gaxpy (parallel gaxpy) 291 - разложение A - xB (?) 420, 421 -- с внешним произведением (outer product) 290 Логарифм матрицы (log of matrix) 490 Локальная программа (nodeprogram) 239 Масштабирование (balancing) 342 Масштабирование по столбцам (column scaling) 30, 121 271 -- строкам (row scaling) 30, 120 Матрица (matrix) - блочная (block) 36 - верхнияя двухдиагональная (upper bidiagonal) 29 -- треугольная (upper triangular) 29 -- хессенбергова (upper Hessenberg) 29 - вырожденная (singular) 29, 30 - диагональная (diagonal) 29, 30 - ленточная (band) 29 - невырожденная (nonsingular) 57 - нижняя двухдиагональная (lower bidiagonal) 29 -- треугольная (lower triangonal) 29 -- хессенбергова (lower bidiagonal) 29 - нормальная (normal) 303 - плохо обусловленная (ill-conditional) 82 - преобразования Гаусса (Gauss transformation) 93 - простая (nonderogatory) 331 - симметричная (symmetric) 33 - строго диагональная доминирующая (diagonal dominance) 116 - треугольная (triangular) 30, 91 -- умножение (multiplication) 92 - трехдиагональная (tridiagonal) 29 - унитреугольная (unit triangolar) 91 - Хаусхолдера (Hausholder) 181, 183 - хорошо обусловленная (well-conditioned) 82 - эрмитова (Hermitian) 36 Матричная норма (matrix norm) 61, 62 -- соглассованность (consistency) 62 -- соотношения (relations between) 64 -- Фробениуса (Frobenius) 62 - экспонента (expnential matrix) 495, 496 -- аппроксимация Паде (Pade approxination) 497 -- чувствительность (sensitivity) 495, 496 Матричные операции (operations with matrices) 17, 18 Машинная точность (machine precision) 67 Метод см.также Алгоритм Метод Аазена (Aasen's method) 152-156 - Арнольди (Arnoldi) 448 - Грамма-Шмидта классический (classical Gram-Schnidt) 202 --- модифицированный (modified) 202 - конечных элементов (finite element method) 284 --- оценка погрешности (error estimation) 293, 294 - Ланцоша (Lanczos) 445-449 - нормальных уравнений (normal estimation) 206 - с диагональным выбором (diagonal pivoting method) 156-158 - сопряженных градиентов (conjugate gradient method) 464-467 --- связь с алгоритмом Ланцоша (?) 468 - Холецкого блочный (block) 136 -- с внешним произведением (outer product) 135 - Холецкого: gaxpy-версия (?) 134 - чебышевских полуитераций (Chebyshev semiatteration method) 457, 458 - Штрассена (Strassen) 42, 43 - Якоби для линейных систем (Jacobi) 453 Методы для неопределенных систем (indefinite system methods) 127 - линейного уравнения для блочных трехдиагональных систем (linear equation method for block tridiagonal systems) 159 ---- ленточных систем (band systems) 141, 142 --- положительно определенных систем (positive definite systems) 132 ---- симметричных неопределенных систем (symmetric indefinite systems) 150 ---- систем Вандермонда (Vandermonde) 166 ---- тёплицевых систем (Toeplitz) () 171 ---- треугольных систем (triangular systems) 86 Минимальная невязка (minimal residual) 206 Многочлен Чебышева (Chebyshev polynomial) 206 Множитель Лагранжа (Lagrange multiplier) 503 Модель с распределенной памятью (destributed memory model) 238, 239 Модификация Краута и Дулитла (Craut and Doolittle) 99 Модификация (update) 21 - QR-разложение (updating QR factorization) 528 Мониторы (monitors) 257, 258 - определенные (?) 257 - cholij 295 - block.chol 296 - nexti 257 - nextij 270 - gax 261 - prod 271 Наибольшее и наименьшее сингулярные числа матрицы (largest and smallest singular matrix numbers) 385 Наискорейший спуск и сопряженные направления (steepest descent and conjugate grandients) 461 Направления спуска (search directions) 462, 463 Невязки метода сопряженных направлений (residuals of conjugate gradient method) 465 - выбор подмножества (subset selection) 512, 513 - задачи наименьших квадратов (LS problem) 512, 513 Недоопределенные системы (underdeterminate systems) 234 Назависимость линейная (linear independence) 56 Ноетрицательно определенные системы (semi-definite systems) 137 Непрерывность разложения собственных значений (cntinuity of eigenvalue dicompositions) 306 Неприводимые матрицы Хессенберга (unreduced Hessenberg matrices) 329 Неравенсто Коши-Шварца (Cauchy-Schwartz inequality) 60 Несимметричная проблема собственных значение (unsymmetric eigenproblem) 299 Нессиметричный метод Ланцоша (unsymmetric Lanczos method) 449 Неявный симметричный QR-шаг со сдвигом Уилкинсона (inplicit symmetric QR step with Wilkinson shift) 380 Норма (norm) 59 Нормальное уравнение (normal equation) 206 Нормальность матрицы и обусловленность собственных значений (normality and eigenvalues condition) 310 Нормальные матрицы (normal matrices) 303 Нормы (norms) - векторные (vector) 59-60 - матричные (matrix) 61 - степеней матрицы (powers of matrix) 321 - Фробениуса (Frobenius) 73-75 Обобщенные сингулярного разложения (generalized singular value decomposition) 502 Обобщенная задача наименьших квадратов с ограничениями (genetalized constrained least squares) 501 Обобщенный метод наименьших квадратов (generalized least squares) 231 Обозначение абсолютной величины (absolute value notation) 67 Обозначения(notation) блочная матрица (block matrix) 37 вектор (vector) 17 двоеточие (colon) 32 матрица (matrix) 17 Образ матрицы (range of a matrix) 57 Обратная матрица (inverse of matrix) 57 -- случай Тёплица (Toeplitz case) 175, 176 - подстановка (back substitution) 88 -- ортоганальная итерация (inverse orthogonal iteration) 322, 323 Обратные задачи о собственных значениях (inverse eigenvalue problems) 532, 533 Обратный метод последовательной верхней релаксации (backward successive over-relaxation) 458 Обращение в машинный нуль (underflow) 66 Обусловленность (condition) 81 - оценка (estimation) 123 - прямоугольных матриц (condition of rectangilar matrix) 206 Одноранговая модификация диагональной матрицы (rank-one modification of a diagonal matrix) 413 -- задачи на собственные значения (?) 542, 545 -- QR-разложений (QR-factorizations) 528 Операции над векторами (vector operations) 18 -- матрицами (matrix operations) 17 Описание алгоритма (specifying algorithm) 27 Ортогональная матрица (orthogonal matrix) 73 Ортогональное дополнение (orthogonal complement) 73 Ортогональные векторы (orthogonal vectors) 73 Ортогональный проектор (orthogonal projection) 77 Ортогональное матричное представление (orthogonal matrix representation) 186 --- вращение Гивенса (Givens rotation) 189 --- факторизованное (factorized form) 185 Ортогональные итерации (orthogonal iterations) 318 Ортонормированный базис (orthonormal basis) 73 Ортонормальность (orthonormality) 73 Основное решение в методе наименьших квадратов (basis solution in LS) 224, 225 Отбрасывание разрядов (cancellation) 66 Отедельность матриц (separation of matrices) 313 Отношение вычислительных затрат к коммуникациям (computions/communications ratio) 246 Оценка прегрешности в степенном методе (error estimation in power method) 317 Ошибка округления (roundoff errors) 68, 70 Параллельные вычисления (parallel computations) -- асинхронная тороидальная процедура (systolic torus) 269 -- матричное умножение, кольцевой алгоритм (matrix multiplication) 273 ---- общая память динамическая (shared memory dynamic) 286 ------ статическая (shared memory static) 286 -- решение трегольной системы (triangelar system solving) 285 ----- кольцевой алгоритм 285 -- Холецкий, асинхронная сеточная процедура (mesh) 285 --- кольцевой алгоритм (?) 273 --- систолический массив (?) 286 -- Якоби циклический метод (cyclic Jacobi) 403 -- gaxpy 238 --- динамическое разделение памяти (?) 259 --- передача сообщений (?) 249 --- статическое разделение памяти (?) 252 --- систолическая модель (?) 243 -- QR-разложение (QR-factorization) 238 --- систолический алгоритм (systolic mesh) 287 Передача сообщений (message passing) 248, 249 Переменные условия (condition variables) 258 Переопределенная система (overdetermined system) 205 Переполнение (overflow) 66 Перестановка циклов (loop reordering) 26, 27 Перестановочные матрицы (permutation matrices) 107 Пересылки в общей памяти (shared memory traffic) 253, 254 Персимметричная матрица (perxymmetric matrix) 172 Поворот пространства (rotations of subspaces) 518 Погрешность абсолютная (absolute error) 60 - матричной функции (matrix function) 488, 489 - относительная (relative) 60 Подматрица (submatrix) 38 Подпространство (subspace) 56 - базис (basis) 56 - инвариантное (invariant) 300, 301 - ортогональная проекция (orthogonal projection) 77 - пересечение (intersection) 522 -- ядер (null-space intersection) 519 - поворот (rotation) 518 - понижающее (deflating) 364 - прямая сумма (direct sum) 56 - расстояние между (distance) 77 - размерность (orthogonal projection) 77 - сингулярное (syngular) 385, 386 - угол между (?) 520, 521 Полиноминальный предобуславливатель (polynomial preconditioner) 477 Полное ортогональное пространство (complete orthogonal space) 217 Полные ортогональные разложения (complete orthogonal decomposition) 217 Положительно определенные системы (positive definite systems) 132, 133, 443 --- алгоритм Ланцоша (Lanczos) 442 --- итерации Гаусса-Зейделя (Gauss-Scidel iterations) 455 --- нессиметричные (unsymmetric) 133 --- свойства (properties) 132, 133 --- симметричные (symmetric) 134 --- LDLt 134 Полярное разложение (polar decomposition) 140 Понижающее подпростраство (deflating subspace) 364 Последовательная верхняя релаксация SOR (successive over-relaxation) 456 --- предобуславливатель (preconditioner) 477 --- симметричная задача собственных значений (symmetric eigenproblem) 364 ---- SSOR 458 --- симметричные неопределенные системы (symmetric endefinite systems) 150 Последовательность Штурма (Sturm sequence) 393 Потеря ортогональности (loss of orthogonality) 437 -- в методе Грама-Шмидта (in Gram-Schmidt) 203 ---- Ланцоша (in Lanczos) 434, 435 - точности (cancellation) 67 Правило Симпсона (Simpson's rule) 494 Предобуславливатели (preconditioners) 473 - неполного разложения Холецкого (incomplete Cholesky) 473 - неполные блочные (incomplete block) 474 Преобразование Гаусса (Gauss transformations) 93 Гаусса-Жордана (Gauss-Jordan) 101-102 - Кэли (Cayley transform) 194 - матриц (transformation of matrices) -- быстрое Гивенса (fast Givens) 191 --, вращение Гивенса (Givens ratations) 191 -- Гаусса (Gauss) 93 -- гиперболические (hyperbolic) 532 -- отражение Хаусхолдера (Hausholder reflections) 182 - подобия (similarity transformation) 300 -- неунитарное (nonunitary) 306 -- определение (definition) 301 -- условия (conditions) 306 Принцип "разделяй и властвуй" ("divide and conquer" principle) 42 Проблема собственных занчений (eigenproblem) 368 --- несимметричная (unsimmetric) 299 --- симметричная (symmetric) 368 Проекция (projections) 77 Простые матрицы (nonderogatory matrices) 332 Процессор id (id processor) 223 Псевдообратный (pseudoinverse) 223 Пучок матриц (pencil) 353 -- диагонализация (diagonalization) 418 -- симметрично-определенный (symmetric-definite) 418 -- собственные значения (eigenvalues) 353 -- эквивалентность (equivalence) 355 Равномерная загруженность (load balancing) 244, 254, 264 Разбиение матрицы по столбцам (partitioning a matrix into columns) 20 --- строкам (rows) 20 Разложение (factorization, decomposition) - жорданово (Jordan decomposition) 305 - Кронекера (Kronecker form) 355 - обобщенное вещественное Шура (generalized real Schur decomposition) 355 -- симметричное 2 х 2 (?) 400 - трехдиагональное (tridiagonal) 376 - Хессенберга (Hessenberga reduction) 326, 327 -- треугольное (Hessenberg triangular) 376 - Холецкого (Cholesky factorization) 134 -- неполное (incomplete) 473 - Шура (Schur decomposition) 302 -- вещественное (?) 325, 326 - QR 195 Размерность (dimension) 56 Ранг матрицы (rank of a matrix) 57 -- выбор подмножества (subset selection) 512 -- численное нахождение (determination) 225 -- численный (numerical) 226 -- QR-разложение (?) 225 -- SVD 75 Распределенные структуры данных (distributed data structures) 242 Расщепление (splitting) 454 Релаксационный параметр (relaxation parameter) 456, 457 Решение задачи наименьших квадратов (least squares solutions) ---- методом быстрых вращение Гивенса (fast Givens) 209, 210 ----- Ланцоша (Lanczos method) 447 ---- модифицированным методом Грамма-Шмидта (modified Gram-Schmidt method) 209 --- преобразования Хаусхолдера (Hausholder reduction) 208 --- SVD 222 - линейной системы (solving a linear system) 97 Свойство чередования (enterlacing property) 369, 370 Сдвиг (shift) - в симметричном QR-алгоритме (?) 380 -- QR-итерации (?) 335 -- QZ-шаге (?) 360 -- SVD-алгоритме (?) 387 - Уилкинсона (Wilkinson) 378 Секулярное уравнение (secular equation) 416, 503 Сетевая топология (network topology) 239 Сети процессоров (processor networks) 239 Сеточная тороидальная процедура (mesh/torus algorithm) --- умножение матриц (matrix times matrix) 269 --- Холецкого алгоритм (?) 285, 286 Сигнатурная матрица (signature matrix) 533 Симметрично-определенный пучок (symmetric-definite pencil) 418 Сингулярное разложение (singular value decomposition (SVD)) 74 -- алгоритм (?) 383, 387 -- доказательство (proof) 74 -- метод Ланцоша (?) 446 --- наименьших квадратов с ограничением (constrained least squares) 503 ----- общая задача (total least squares) 514 -- обобщенное (generalized) 422 -- пересечение подпростраств (subspace intersection) 522 -- поворот подпространств (subspace rotation) 518, 519 -- проекции (projection) 77 -- псевдообратная матрица (pseudo inverse) 223 -- ранг матрицы (rank of a matrix) 75 -- ядро (null-space) 75 Сингулярные значения (singular values) 74 -- возмущения (perturvations) 384 -- матрицы (singular values of a matrix) 74 -- минимаксная характеризация (minimax characterization) 369, 370 -- собственные значения (eigenvalues) 383 Сингулярный вектор левый (?) 74 -- нуль-пространство (null-space) 75 -- образ (range) 75 -- правый (?) 74 Синус матрицы (sine of a matrix) 490 Системы Вандермонда (Vandermonde systems) 166, 167 - общей памяти (shared memory systems) 253, 254 Скалярное произведение (dot product) 18 След (trace) 300 Сложение (addition) 17 Собственное значение -- внутреннее (interior) 432 -- дефектное (defective eigenvalue) 305 Собственные значения (eigenvalues) 305 -- детерминант (determinant) 317 -- доминирующее (dominant) 317 -- комплексно сопряженные (complex conjugate) 325, 326 -- обобщенные (genetalized) 354 -- последовательность Штурма (Sturm sequence) 369 -- симметрические матрицы (symmetric matrices) 369 -- след (trace) 300 -- упорядочение в вещественной форме Шура (ordering in Schur form) -- характеристический многочлен (characteristic polynomial) 299, 300 Собственный вектор (eigenvector) 300 -- возмущение (perturbation) 311 -- дефектная матрица (defective matrix) 306 -- левый (left) 300 -- плохо обусловленная матрица (ill-conditioned matrix) 306 -- правый (right) 300 -- формула Шура (Shur form) 403 -- чувствительность (sensitivity) 311 Сопряженное транспонирование (conjugate transposition) 27 Сопряженные направления (conjugate directions) 463 Сосед (neighbor) 240 Спектр (spectrum) 300 Спектральный радиус (spectral radius) 454 Стационарные значения (stationary values) 369 -- с ограничениями (constrained) 523 Степени матрицы (matrix powers) 381 Сепенной метод (power method) 316 - ряд матрицы (power series of a matrix) 490 Стоимость обмена информацией (communication cost) 246 Структура данных для блочных матриц (block data structure) 41, 42 -- ленточной матрицы (vector computers) 143 -- при векторной обработке (vector computers) 51, 52 Сходимость (convergence) - итераций Якоби (Jacobi iterations) 344 - итерационных методов (iterative methods) 454 - метода бисекций (bisection) 393 -- Гаусса-Зейделя (Gauss-Seidel) 455 -- итераций с отношениями Рэлэя (Rayleigh quotient iterations ) () 396 -- Ланцоша (Lanczos) 430 -- наискорейшего спуска (steepest descent) 461 -- сопряженных градиентов (conjugate gradients) 468, 469 -- Якоби для симметричной задачи собственных значений (Jcobi's method for the symmetric eigenproblem) 401 - чебышевских полуитераций (Chebyshev semi-iterative) 457 - обратных итераций (inverse iterations) 344 - ортогональной итерации (orthogonal iteration) 318 - симметричного QR-шага со сдвигом (symmetric QR iteration) 380 - степенных итераций (power method) 316 - циклического метода Якоби (cyclic Jocobi) 403 - QR-алгоритма (QR) 335 - QZ-алгоритма (QZ) 353 - SVD-алгоритма (SVD) 391 Теорема (theorem) Бауэра-Файка (Bauer-Fike) 308 Виландта-Хоффмана о собственных значениях (Wielandt-Hoffman) 370 Куранта-Фишера о минимаксе (Courant-Fischer minimax theorem) 369 о кругах Гершгорина (Gershgorin circle theorem) 307, 308 - минмаксе для собственных значений (?) 369 - неявном Q (implicit Q theorem) 330, 331 Сильвестра об инерции (Sylvester law of inertia) 374 сингулярных значений (singular values) 385 Теории возмущения (pertubation theory) 307 -- для инвариантных подпространств (invariant subspaces) 313 -- для недоопределенных систем (underdetermined systems) 236 --- псевдообратной матрицы (pseudo inverse) 223 --- собственных векторов (eigenvectors) 310, 311 ---- значений (eigenvalues) 307 -- инвариантные подпространства симметричных матриц (invariant subspaces of symmetric matrices) 371 -- линейной системы (linear equation problem) 80, 81 -- обобщенное собственное значение (generalized eigenvalues) 355 -- пара сингулярных подпространств (sigular subspace paire) 385, 386 -- сингулярные значения (singular values) 384, 385 -- собственные значения симметричной матрицы (eigenvalue of symmetric matrix) 370 - Каниеля-Пейджа (Kaniel-Paige threory) 430 Тёплицева матрица (Toeplitz matrix) 171 - система (Toeplitz system) 171 Тор (torus) 240 Точность (precision) 66 Треугольные системы (triangular systems) -- ленточные (band) 142, 143 -- неквадратные (non-square) 91 -- случай нескольких правых частей (multiple) 89 Трехдиагональная матрица (tridiagonal matrix) 475 -- обратная (inverse) 475 Трехдиагональные системы (tridiagonal systems) 146 Трехдиагонализация (tridiagonalization) 376, 427 - Ланцоша (Lanczos) 427, 428 -- блочный вариант (block version) 438 -- внутренние собственные значения (interior eigenvalues) 432 -- с выборочной ортогонализацией (selective orthogonalization) 437 --- полной переортогонализацией (complete reothogonalization) 436 -- s-шаговый (s-step) 440 -- сопряженные градиенты (conjugate gradients) 439 -- степенной метод (power method) 431 - подпространство Крылова (Krylov subspace) 426 - Хаусхолдера (Hausholder) 377 Умножение блочных матриц (block matrix multiplication) 39 - матриц по Штрассену (Strassen multiplication) 43 ---- блочная версия (block version) 19 - матрицы на вектор (matrix-vector) 20 --- матрицу (matrix-matrix) 17, 23 --- число (scalar matrix multiplication) 39 Уравнение Сильвестра (Sylvester equation) 347 Уравновешивание строчностолбцевое (row-column equilibration) 121 Уровень (level) 23, 27 - операции (of operation) 23 Ускорение (speed-up) 246 Ускорение Ритца (Ritz acceleration) 396, 397 -- пара и метод Ланцоша (pair and Lanczos method) 437, 438 Условия Мура-Пенроуза (Moore-Penrose condition) 223 Устоцчивость (stability) 498, 499 Флоп (flop) 30, 31 Формула двойных углов (doubling formulae) 492, 497 - Шермана-Моррисона (Sherman-Morrison formula) 57 Функции от матриц (matrix functions) 482 --- интегральные (integrating) 494 --- многочлен (polynomial equation) 482, 493 Функция треугольной матрицы (function of triangelar matrix) 484, 485, 488 Характеристический многочлен (characteristic polynomial) 299 -- обобщенная задача собственных значений (generalized eigenproblem) 354 Хессенбергова форма и итерация Арнольди (Arnoldi process and Hessenberg form) 488 --- преобразования Хаусхолдера (Haushelder reduction) 328 --- QR-итерация (?) 344 --- QR-разбиения (?) 199 -- неприводимость матрицы (unreduced) 330 -- обратные итерации (inverse) 328 -- свойства (?) 329 Хессенбергово-треугольное преобразование (Hessenberg triangular form reduction) 357, 358 Хессенберговы системы (Hessenberg systems) 145, 146 Хранение ленточных матриц (band matrix store) 32, 33 - по диагоналям () 33, 34 Циклическая редукция (cyclic reduction) 162 Циклический метод Якоби (cyclic Jacobi method) 402, 403 Число Кроуфорда (Crowford number) 420 - обусловленности (condition number) 81, 82 -- по Шкеелю (Skeel condition number) 85 - с плавающей точкой (floating point number) 65 Чистка (sweep) 402 Чувствительность. См. Теория возмущения Чувствительность линейного уравнения (linear equation sensitivity) 83-85 Ширина ленты (bandwidth) 30 -- верхняя (upper bandwidth) 30 -- нижняя (lower bandwidth) 30 Штрассена алгоритм (Strassen algorithm) 42 Эквивалентность норм (equivalence of norms) 63 Эффективность (efficiency) 247 Ядро матрицы (null-space) 57 - пересечение (intersection) 519 A-норма (A-norm) 526 A-сопряженное направление спуска (A conjugate searching direction) 463 cols 20 continue 259 CS-разложение (CS decomposition) 79 delay 259 gauss 93 gaxpy 22 - алгоритм (gaxpy algorithms) gaxpy в распределенной памяти (in distributed memory) 243, 246, 248, 250 - на общей памяти (in shared memory) 252, 255, 258, 261 - и внешние произведения (gaxpy vs. outer product) 52 get 254 givens 188 glob. init 253 house 183 LDLt-разложение 130, 131 LDMt-разложение 127, 128 length 18 loc.init 241 LR-итерация (LR iterations) 321 LU-разложение (LU factorization) - ленточное (band) 141, 142 - дифференцирование (differentiation) 101 - детерминанта (determinant) 95 - существование (existence) 96, 97 Mathlab 18 mesh 239 my. id 241 n x n матрица перестановок (n x n permitation matrix) 172 p-нормы (p-norms) 60 - минизация (minimisation in) 205 put 254 quit 241 QR-алгоритмы для собственных значений (QR-iterations) 316 QR-итерации (QR-iterations) - вывод (derivation) 319, 320 - неявный двойной сдвиг (implicit double shift) 338 - неявный сдвиг (?) 337 - одинарный сдвиг (single shift) 336 - сдвиг и точное собственное значение (?) 336 QR-разложение (QR-decomposition) 195 - вычисление методом блочных отражений (block Hausholder computation) 197 -- классичеким методом Грама-Шмидта (classical Gram-Schmidt) 201 -- методом быстрых вращений Гивенса (fast Givens rotations) 200 -- модифицированным методом Грамма-Шмидта (modified) 202 - квадратные системы (square systems) 235 - модификация (updating) 528 - наименьших квадратов (LS) 208 - неполного ранга с выбором ведущего столбца (rank defficient column pivoting) 215 - преобразования Хаусхолдера (Hausholder computation) 196, 197 - ранг матрицы (rank of matrix) 215, 216 - свойства (?) 201 - сеточная систолическая процедура (systolic mesh) 287 - симметричное (symmetric) 376 - хессенберговой матрицы 199 QR-шаг Фрэнсиса (Francis QR-step) 340 R-двухдиагонализация (R-bidiagonalization) 219 recv 214 saxpy 18 send 241 span 57 SVD-шаг Голуба-Кахана (Golub-Kahan SVD step) 389 sym.schur 401 WY-представление (WY-representation) 187 2-норма Глава 1. Умножение матриц 16 1.1 Основные алгоритмы и обозначения 16 1.1.1. Обозначения для матриц 1.1.2. Операции над матрицами 1.1.3. Обозначения для векторов 1.1.4. Операции над векторами 1.1.5. Вычисление скалярных произведений и операций saxpy 1.1.6. О формализме и стиле изложения 1.1.7. Умножение матрицы на вектор 1.1.8. Разбиение матрицы на строки и столбцы 1.1.9. Использование дваеточия 1.1.10. Мидификация внешним произведением 1.1.11. Вычисление операции gaxpy 1.1.12. Понятие уровня 1.1.13. Умножение матриц 1.1.14. Умножение матриц с использованием скалярных произведений 1.1.15. Умножение матриц с использованием gaxpy 1.1.16. Умножение матриц с использованием внешних произведений 1.1.17. Перестановка циклов 1.1.18. О матричных соотношениях 1.1.19. Об описании алгоритмов 1.1.20. Комплексные матрицы Задачи. Замечания и литература. 1.2 Учет структуры матрицы 29 1.2.1. Ленточные матрицы и обозначения для них 1.2.2. Действия с диагональными матрицами 1.2.3. Умножение треугольных матриц 1.2.4. Флопы 1.2.5. Еще раз об использовании дваеточия 1.2.6. Хранение ленточных матриц 1.2.7. Симметрии 1.2.8. Хранение по диагоналям 1.2.9. Запись поверх входных данных 1.3 Блочные матрицы и алгоритмы 36 1.3.1. Обозначения для блочных матриц 1.3.2. Операции с блочными матрицами 1.3.3. Обозначения для подматриц 1.3.4. Умножение блочной матрицы на вектор 1.3.5. Перемножение блочных матриц 1.3.6. Важный частный случай 1.3.7. Структуры данных для блочных матриц 1.3.8. Умножение матриц, основанное на принципе "разделяй и властвуй" 1.4 Некоторые аспекты векторно-ковейерных вычислений 45 1.4.1. Конвейеризация арифметических операций 1.4.2. Векторные операции 1.4.3. Длина вектора 1.4.4. Оптимизация кода с учетом длины векторов 1.4.5. Многоуровневая память 1.4.6. Единичный шаг выборки 1.4.7. Роль структур данных 1.4.8. Операции , внешние произведения и вектоные обмены 1.4.9. Блочные алгоритмы, кэш-память и повторное использование данных 1.4.10. Резюме Глава 2. Матричный анализ 56 2.1 Основные сведения из линейной алгебры 56 2.1.1. Линейная независимость, подпространства, базис и размерность 2.1.2. Область значений, ядро и ранг матрицы 2.1.3. Обратная матрицы 2.1.4. Детерминант 2.1.5. Дифференцироание 2.2 Векторыне нормы 59 2.2.1. Определения 2.2.2. Некоторые свойства векторных норм 2.2.3. Абсолютная и относительная погрешности 2.2.4. Сходимость 2.3 Матричные нормы 61 2.3.1. Определения 2.3.2. Некоторые свойства матричных норм 2.3.3. Матричная 2-норма 2.3.4. Возмущения и обратная матрица 2.4 Матричные вычилсления с конечной точностью 65 2.4.1. Числа с плавающей точкой 2.4.2. Модель арифметики с плавающей точкой 2.4.3. Потеря точности 2.4.4. Обозначение абсолютный величины 2.4.5. Ошибки округления в скалярных произведениях 2.4.6. Альтернативные методы количественной оценки ошибок 2.4.7. Вычисление скалярных произведений с накоплением 2.4.8. Ошибки округления в других основных матричных вычислениях 2.4.9. Прямой и обратный анализ ошибок 2.4.10. Ошибки округления в алгоритме Штрассена 2.5 Ортогональность и сингулярное разложение 73 2.5.1. Ортогональность 2.5.2. Нормы и ортогональные преобразования 2.5.3. Сингулярное разложение 74 2.5.4. Неполнота ранга и SVD 2.6 Проекции и CS-разложение 77 2.6.1. Ортогональные проекторы 2.6.2. Проекторы, связанные с SVD 2.6.3. Расстояния, связанные в подпространствами 2.7 Чувствительность квадратных систем к возмущениям 80 2.7.1. SVD-анализ 2.7.2. Обусловленность 2.7.3. Определители и близость к вырожденности 2.7.4. Точная оценка по норме 2.7.5. Несколько точных покомпонентных оценок Глава 3. Линейные системы общего вида 87 3.1 Треугольные системы 87 3.1.1. Прямая подстановка 3.1.2. Обратная подстановка 3.1.3. Столбцовые версии 3.1.4. Случай нескольких правых частей 3.1.5. Доля флопов 3 уровня 3.1.6. Решение неквадратных треугольных систем 3.1.7. Унитреугольные системы 3.1.8. Алгебраические свойства треугольных матриц 3.2 LU-разложение 92 3.2.1. Матрица преобразования Гаусса 3.2.2. Применение матриц преобразования Гаусса 3.2.3. Свойства ошибок округления в преобразованиях Гаусса 3.2.4. Приседение к верхнему треугольному виду 3.2.5. LU-разложение 3.2.6. Несколько практических замечаний 3.2.7. Где хранить матрицк L? 3.2.8. Решение линейной системы 3.2.9. Gaxpy-верси LU-разложения 3.2.10. Модификация Краута-Дулитла 3.2.11. Блочное LU-разложение 3.2.12. LU-разложение прямоугольной матрицы 3.2.13. Несостоятельность метода 3.3 Анализ ошибок округления в методе исключения Гаусса 102 3.3.1. Ошибки в LU-разложении 3.3.2. Решение треугольных систем с приближенными треугольными матрицами 3.4 Выбор ведущего элемента 106 3.4.1. Перестановочные матрицы 3.4.2. Частичный вабор ведущего элемента: общая идея 3.4.3. Детали стратегии частичного выбора 3.4.4. Где хранить матрицу L? 3.4.5. Gaxpy-версия 3.4.6. Анализ ошибок округления 3.4.7. Метод блочного исключения Гаусса 3.4.8. Полный выбор ведущего элемента 3.4.9. Замечания по стратегии с полным выбором ведущего элемента 3.4.10. Отказ от выбора ведущего элемента 3.4.11. Некоторые приложения 3.5 Уточнение и оценивание точности 119 3.5.1. Большая невяжка дает плохую точность 3.5.2. Масштабирование 3.5.3. Итерационное уточнение 3.5.4. Оценка обусловленности Глава 4. Линейные системы специального вида 127 4.1 Разложение вида LDM LDL 127 4.1.1. LDMt-разложение 4.1.2. Симметрия и LDLt-разложение 4.2 Положительно определенные системы 132 4.2.1. Положительная определенность 4.2.2. Несимметричные положительно определенные системы 4.2.3. Симметричные полжительно определенные системы 4.2.4. Gaxpy-версия разложения Холецкого 4.2.5. Метод Холецкмого с внешним произведением 4.2.6. Блочный метод Холецкого со скалярным произведением 4.2.7. Устрйчивость процесса Холецкого 4.2.8. Неотрицательно определенные матрицы 4.2.9. Симметричный выбор ведущего элемента 4.3 Ленточные системы 141 4.3.1. Ленточное LU-разложение 4.3.2. Решение треугольных ленточных систем 4.3.3. Структура данных летночной матрицы 4.3.4. Ленточное исключение Гаусса с выбором ведущего элемента 4.3.5. LU-разложение матрицы Хессенберга 4.3.6. Ленточный метод Холецкого 4.3.7. Решение трехдиагональных систем 4.3.8. Вопросы векторизации 4.4 Симметричные неопределенные системы 150 4.4.1. Алгоритм Парлетта-Рейда 4.4.2. Метод Аазена 4.4.3. Выбор ведущего элеметна в методе Аазена 4.4.4. Методы с диагональным выбором 4.4.5. Устойчивость и эффективность 4.4.6. Сравнение метода Аазена с диагональным выбором 4.5 Блочные трехдиагональные системы 159 4.5.1. Блочное LU-разложение 4.5.2. Блочное диагональное доминирование 4.5.3. Сравнение блочного и ленточного решений 4.5.4. Блочная циклическая редукция 4.6 Системы Вандермонда 166 4.6.1. Интерполяционный многочлен V(T)a = f 4.6.2. Система Vz = b 4.6.3. Устойчивость 4.7 Теплицевы системы 171 4.7.1. Три задачи 4.7.2. Решение уравнений Юла-Уолкера 4.7.3. Задача с произвольной правой частью 4.7.4. Вычисление обратной матрицы 4.7.5. Вопросы устойчивости Глава 5. Ортогонализация и метод наименьших квадратов 181 5.1 Матрицы Хаусхолдера и Гивенса 181 5.1.1. Двумерный случай 5.1.2. Отражение Хаусхолдера 5.1.3. Вычисление вектора Хаусхолдера 5.1.4. Умножение на матрицы Хаусхолдера 5.1.5. Ошибки округления 5.1.6. Факторизованное представление 5.1.7. Блочное представление 5.1.8. Вращение Гивенса 5.1.9. Умножения на матрицы Гивенса 5.1.10. Ошибки округления 5.1.11. Представление произведений матриц Гивенса 5.1.12. Распространение ошибок 5.1.13. Быстрые вращения 5.2 QR-разложение 195 5.2.1. QR-разложение: преобразования Хаусхолдера 5.2.2. QR-разложение: метод блочных отражений 5.2.3. QR-разложение: преобразования Гивенса 5.2.4. QR-разложение хессенберговской матрицы с использоанием вращений 5.2.5. QR-разложение: быстрые вращения 5.2.6. Свойства QR-разложения 5.2.7. Классический метод Грама-Шмидта 5.2.8. Модифицированный метод Грама-Шмидта 5.2.9. Вычислительные затраты и точность 5.3 Задача наименьших квадратов: случай полного ранга 205 5.3.1. Следствие полноты ранга 5.3.2. Обусловленность прямоугольных матриц 5.3.3. Метод нормальных уравнений 5.3.4. Решение задач LS через QR-разложение 5.3.5. Срыв в случае "почти неполноранговости" 5.3.6. О методе MGS 5.3.7. Решение задачи LS методом быстрых вращений 5.3.8. Чувствиетльность задачи LS к возмущениям 5.3.9. Сравнение метода нормальных уравнений и QR-разложения 5.4 Другие ортогональные разложения 215 5.4.1. Случай неполного ранга: QR с выбором ведущего столбца 5.4.2. Полные ортогональные разложения 5.4.3. Двухдиагонализация 5.4.4. R-двухдиагонализация 5.4.5. Сингулярное разложение и его вычисление 5.5 Задача LS неполного ранга 221 5.5.1. Решение с минимальной нормой 5.5.2. Полное ортогональное разложение и x(ls) 5.5.3. Сингулярное разложение и задача LS 5.5.4. Псевдообратная матрица 5.5.5. Некоторые вопросы чувствительности к возмущениям 5.5.6. QR с выбором ведущего столбца и основные решения 5.5.7. Численное нахождение ранга с помощью aП = QR 5.5.8. Численный ранг и SVD 5.5.9. Некоторые сравнения 5.6 Взвешивание и итерационное уточнение 230 5.6.1. Взвешивание по столбцам 5.6.2. Взвешивание по строкам 5.6.3. Обобщенные наименьшие квадраты 5.6.4. Итерационное уточнение 5.7 Квадратные и недоопределенные системы 234 5.7.1. Использование QR и SVD систем для решения квадратных систем 5.7.2. Недоопределенные сисмемы 5.7.3. Возмущенные недоопределенные системы Глава 6. Параллельные матричные вычисления 6.1 Операции на распределенной памяти 238 6.1.1. Системы с распределенной памятью 6.1.2. Сети процессоров 6.1.3. Имена для соседей 6.1.4. Инициализация и окончание 6.1.5. Связь между процессорами 6.1.6. Некоторые распределенные структуры данных 6.1.7. Систолическая модель 6.1.8. Атрибуты параллельного алгоритма 6.1.9. Равномерная загруженность 6.1.10. Стоимость обмена информацией 6.1.11. Эффективность и ускорение 6.1.12. Дробление вычислений 6.1.13. Модель передачи сообщений 6.1.14. Дальнейшие примеры передачи сообщений 6.2 Операции на общей памяти 252 6.2.1. Статическое планирование операции gaxpy 6.2.2. Обмен данными с общей памяти 6.2.3. Равномерная загруженность 6.2.4. Синхронизация с помощью барьеров 6.2.5. Парадигма динамического резерва задач 6.2.6. Мониторы 6.2.7. Столбцовая реализация операции gaxpy 6.3 Параллельное умножение матриц 262 6.3.1. Процедуры для блочной операции 6.3.2. Проблема равномерной загруженности 6.3.3. Некоторые вопросы, связанные с дроблением вычислений 6.3.4. Систолическое умножение 3 x 3-матриц 6.3.5. Общий случай 6.3.6. Блочный аналог 6.3.7. Асинхронная тороидальная процедура 6.3.8. Использование блочного скалярного произведения в режиме резерва задач 6.4 Кольцевые процедуры разложения 273 6.4.1. Кольцевой алгоритм Холецкого: n = p 6.4.2. Кольцевой метод Холецкого: общий случай 6.4.3. Кольцевые процедуры для других разложений 6.4.4. Параллельное решение треугольных систем 6.5 Сеточные процедуры разложения 281 6.5.1. Сеточная систолическая процедура Холецкого 6.5.2. Асинхронная процедура Холецкого 6.5.3. Сеточная систолическая процедура QR-разложения 6.6 Методы разложения на общей памяти 290 6.6.1. Статическое планирование для Холецского с внешними произведениями 6.6.2. Статическое планирование для других разложений 6.6.3. Две параллельные реализации для Холецкого с gaxpy 6.6.4. Замечания об обменах с общей памятью 6.6.5. Блочный Холецкий на базе резерва задач Глава 7. Несимметричная проблема собственных значений 299 7.1 Свойства и разложения 300 7.1.1. Собственные значения и инвариантные подпространства 7.1.2. Основные унитарные разложения 7.1.3. Неунитарные преобразования 7.1.4. Некоторые замечания о неунитарном подобии 7.2 Теория возмущения 307 7.2.1. Чувствительность собственного значения 7.2.2. Обусловленность простого собственного значения 7.2.3. Чувствительность кратных собственных значений 7.2.4. Чувствительность собственного вектора 7.2.5. Чувствительность инвариантного подпространства 7.3 Степенные итерации 316 7.3.1. Степенной метод 7.3.2. Ортогональные итерации 7.3.3. QR-итерации 7.3.4. LR-итерации Приложение 7.4 Хессенбергова форма и вещественная форма Шура 324 7.4.1. Вещественное разложение Шура 7.4.2. Хессенбергов QR-шаг 7.4.3. Разложение Хессетберга 7.4.4. Трехуровневые варианты 7.4.5. Важные свойства матрицы Хессенберга 7.4.6. Сопровождающая матрица 7.4.7. Хессенбергово преобразование через преобразование Гаусса 7.5 Практический QR-алгоритм 334 7.5.1. Исчерпывание 7.5.2. QR-итерации со сдвигами 7.5.3. Стратегия одинарного сдвига 7.5.4. Стратегия двойного сдвига 7.5.5. Стратегия неявного двойного сдвига 7.5.6. Полный процесс 7.5.7. Масштабирование 7.6 Методы вычисления инвариантных подпространств 343 7.6.1. Выделение собственных векторов при помощи обратных итераций 7.6.2. Упорядочение собственных значений в вещественной форме Шура 7.6.3. Блочная диагонализация 7.6.4. Базис собственных векторов 7.6.5. Определение блочных жордановых структур 7.7 QZ-метод для Ax = lyambdaBx 353 7.7.1. Основы теории 7.7.2. Обобщенное разложение Шура 7.7.3. Выводы о чувствительности 7.7.4. Хессенбергово-треугольная форма 7.7.5. Исчерпывание 7.7.6. QZ-шаг 7.7.7. QZ-процесс в целом 7.7.8. Вычисление обобщенных инвариантных подпространств Глава 8. Сииметрична проблема проблема собственных значений 368 8.1 Математические основы 368 8.1.1. Собственные значения симметричных матриц 8.1.2. Теорема о минимаксе и некоторые следствия 8.1.3. Дополнительные результаты 8.1.4. Чувствительность инвариантных подпространств 8.1.5. Закон инерции 8.2 Симметричный QR-алгоритм 376 8.2.1. Преобразование к трехдиагональному виду 8.2.2. QR-итерации с явным одинарным сдвигом 8.2.3. Вариант с неявным сдвигом 8.3 Вычисление SVD 383 8.3.1. Теория возмущения и свойства 8.3.2. SVD-алгоритм 8.4 Некоторые специальные методы 392 8.4.1. Бисекция 8.4.2. Итерация с отношением Рэлея 8.4.3. Ортогональные итерации с ускорением Ритца 8.5 Методы Якоби 399 8.5.1. Идея метода Якоби 8.5.2. Симметричное разложение Шура размера 2 x 2 8.5.3. Обновления в методе Якоби 8.5.4. Классический алгоритм Якоби 8.5.5. Алгоритм циклический по строкам 8.5.6. Барьерный метод Якоби 8.5.7. Анализ ошибок округления 8.5.8. Сравнение с симметричным QR-алгоритмом 8.5.9. Параллельное упорядочение 8.5.10. Кольцевая процедура 8.5.11. Блочная процедура Якоби 8.5.12. SVD-процедура Якоби 8.6 Метод разделяй и властвуй 412 8.6.1. Расщепление 8.6.2. Объединение разложений Шура 8.6.3. Собственная система матрицы D + pzz(t) 8.6.4. Практический синтез 8.6.5. Полный процесс с параллелизмом 8.7 Более общие проблемы собственных значений 418 8.7.1. Математические основы 8.7.2. Методы для симметрично-определенной проблемы 8.7.3. Обобщенная проблема сингулярных значений Глава 9. Методы Ланцоша 426 9.1 Выводы свойства сходимости 426 9.1.1. Подпространства Крылова 9.1.2. Трехдиагонализация 9.1.3. Окончание и оценки погрешностей 9.1.4. Теория сходимости Каниэля-Пейджа 9.1.5. Сравнение степенного метода с методом Ланцоша 9.1.6. Сходимость внутренних собственных значений 9.2 Практические процедуры Ланцоша 433 9.2.1. Реализация в точной арифметике 9.2.2. Ошибки округления 9.2.3. Метод Ланцоша с полной переортогонализацией 9.2.4. Выборочная ортогонализация 9.2.5. Проблема теневых собственных значений 9.2.6. Блочный Ланцош 9.2.7. s-Шаговый алгоритм Ланцоша 9.3 Приложения и обобщения 442 9.3.1. Симметричные положительно определенные системы 9.3.2. Симметричные неопределенные системы 9.3.3. Двухдиагонализация и SVD 9.3.4. Наименьшие квадраты 9.3.5. Идея Арнольди 9.3.6. Несимметричная трехдиагонализация Ланцоша Глава 10. Итерационные методы для линейных систем 452 10.1 Стандартные итерации 452 10.1.1. Итерации Якоби и Гаусса-Зейделя 10.1.2. Расщепление и сходимость 10.1.3. Практическая реализация Гаусса-Зейделя 10.1.4. Последовательная верхняя релаксация 10.1.5. Метод Чебышевских полуитераций 10.1.6. Симметричный метод SOR 10.2 Методы сопряженных градиентов 461 10.2.1. Наискорейший спуск 10.2.2. Произвольные направления спуска 10.2.3. A-сопряженные направления спуска 10.2.4. Метод сопряженных градиентов 10.2.5. Несколько крайне необходимых наблюдений 10.2.6. Связь с алгоритмом Ланцоша 10.2.7. Некоторые практические детали 10.2.8. Сходимость метода 10.3 Сопряженные градиенты с предобусловливанием 471 10.3.1. Вывод 10.3.2. Предобусловливатели, связанные с неполным разложением Холецского 10.3.3. Неполные блочные предобусловливатели 10.3.4. Идеи декомпозиции области 10.3.5. Полиноминальные предобусловливатели 10.3.6. Заключительное предложение Глава 11. Функции от матриц 482 11.1 Спектральные методы 482 11.1.1. Определение 11.1.2. Жорданова характеризация 11.1.3. Метод разложения Шура 11.1.4. Метод блочного разложения Шура 11.2 Аппроксимационные методы 488 11.2.1. Жорданов анализ 11.2.2. Анализ на основе разложения Шура 11.2.3. Тейлоровы апроксиманты 11.2.4. Оценка многочлена от матрицы 11.2.5. Вычисление степеней матрицы 11.2.6. Интегральные функции от матрицы 11.3 Матричная экспонента 495 11.3.1. Теория возмущений 11.3.2. Метод аппроксимации Паде 11.3.3. Некоторые выводы об устойчивости Глава 12. Специальные разделы 501 12.1 Задача наименьших квадратов с ограничениями 501 12.1.1. Задача LSQI 12.1.2. LS-минимизация на сфере 12.1.3. Гребневая регрессия 12.1.4. Задача наименьших квадратов с ограничениями типа равенств 12.1.5. Метод взвешивания 12.2 Выбор подмножеств при помощи SVD 509 12.2.1. Метод QR со столбцовым выбором 12.2.2. Использование метода SVD 12.2.3. Еще раз о противоречии между независимостью столбцов и невязкой 12.3 Общая задача наименьших квадратов 514 12.3.1. Математическая основа 12.3.2. Вычисления для случая k = 1 12.3.3. Геометрическая интерпретация 12.4 Сравнение подпространств при помощи SVD 518 12.4.1. Поворот подпространств 12.4.2. Пересечение ядер 12.4.3. Углы между подпространствами 12.4.4. Пересечение подпространств 12.5 Модифицированные задачи на собственный значения 523 12.5.1. Стационарные значения квадратичной формы с ограничениями 12.5.2. Обратная задача на собственные значения 12.5.3. Одноранговая модификация задачи на собственые значения 12.5.4. Вторая обратная задача на собственные значения 12.6 Модификация QR-разложения 528 12.6.1. Одноранговые модификации 12.6.2. Добавление или удаление столбца 12.6.3. Добавление или исключение строки 12.6.4. Методы гиперболического преобразования